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ABSTRACT 

A  numerical  method  is  presented  for  analyzing  the  transonic 
potential  flow  past  a  lifting,  swept  wing.   A  finite-difference 
approximation  to  the  full  potential  equation  is  solved  in  a  coordi- 
nate system  which  is  nearly  conformally  mapped  from  the  physical 
space  in  planes  parallel  to  the  symmetry  plane,  and  reduces  the 
wing  surface  to  a  portion  of  one  boundary  of  the  computational  grid. 
A  coordinate  invariant,  rotated  difference  scheme  is  used,  and  the 
difference  equations  are  solved  by  relaxation.   The  method  is 
capable  of  treating  wings  of  arbitrary  planform   and  dihedral, 
although  approximations  in  treating  the  tips  and  vortex  sheet  make 
its  accuracy  suspect  for  wings  of  small  aspect  ratio.   Comparisons 
of  calculated  results  with  experimental  data  are  shown  for  examples 
of  both  conventional  and  supercritical  transport  wings.  Agreement 
is  quite  good  for  both  types,  but  it  was  found  necessary  to  account 
for  the  displacement  effect  of  the  boundary  layer  for  the  super- 
critical wing,  presumably  because  of  its   greater  sensitivity  to 
changes  in  effective  geometry. 


INTRODUCTION 

The  development  of  profile  shapes  capable  of  efficient 
operation  in  the  transonic  regime  has  spurred  interest  in  flight 
vehicles  designed  specifically  to  operate  at  near  sonic  speeds. 
The  ability  to  predict  accurately  the  aerodynamic  characteristics 
of  the  complete  three-dimensional  wing  should  have  a  substantial 
impact  on  the  design  of  such  vehicles  by  allowing  detailed  trade- 
off studies  to  be  performed  without  recourse  to  wind  tunnel  test- 
ing of  every  design  variation. 

Recent  advances  in  the  theoretical  prediction  of  inviscid 
transonic  flow  fields  are  based  largely  on  type-dependent,  finite- 
difference  solutions  of  the  steady  potential  equation.  These 
methods  were  first  applied  to  the  transonic  small  disturbance 
equation  by  Murman  and  Cole  [1],  and  the  full  potential  equation 
by  Jameson  [2]  and  Garabedian  and  Korn  [3]  for  the  prediction  of 
airfoil  flow  fields.   The  three-dimensional  small  disturbance 
equation  has  also  been  solved  for  swept  wings  by  Ballhaus  and 
Bailey  [4]  and  for  wing -cylinder  combinations  by  Bailey  and 
Ballhaus  [5].   Finally,  the  full  potential  equation  has  been 
solved  by  Jameson  for  the  transonic  flow  over  an  oblique  yawed 
wing  [6].   Although  an  oblique  wing  should  be  aerodynamically  more 
efficient  than  a  conventional   swept  wing  [7] ,  it  presents 
problems  of  stability  and  control  and  aeroelastic  divergence. 
We  consider  here  the  prediction  of  the  flow  over  a  swept  wing. 

In  Jameson's  treatment  of  the  flow  over  oblique  wings, 
the  coordinate  system  is  aligned  in  planes  normal  to  the  wing 


leading  edge.   Thus,  for   nonzero  angles  of  yaw  the  free  stream 
velocity  vector  is  not  contained  in  these  planes,  and  the  treat- 
ment of  a  symmetry  plane  in  the  flow  past  a  swept  wing  would  be 
difficult  in  this  coordinate  system.   In  the  analysis  presented 
here,  the  flow  is  analyzed  in  coordinate  planes  parallel  to  the 
free  stream  velocity  vector,  and  the  symmetry  condition  is 
applied  on  a  single  coordinate  surface.   To  allow  the  use  of  a 
fine  mesh  to  resolve  the  details  of  the  flow  in  the  sensitive 
region  near  the  leading  edge,  the  spanwise  coordinate  lines  are 
aligned  with  the  leading  edge.   Thus  for  wings  of  appreciable 
sweep,  the  resulting  coordinate  system  is  highly  nonorthogonal . 

The  type  of  geometry  we  shall  treat  is  illustrated  in 
Figure  1.   It  consists  of  a  wing  of  arbitrary  planform   and 
dihedral  extending  from  a  symmetry  plane  (or  wall) .   We  shall 
solve  a  finite  difference  approximation  to  the  full  potential 
equation  for  the  transonic  flow  past  such  a  configuration  using 
a  generalized  relaxation  method.   The  finite  difference  approxi- 
mation is  the  rotated  difference  scheme  introduced  by  Jameson  [6] , 
and  is  not  in  conservation  form.   This  can  introduce  substantial 
errors  in  the  treatment  of  flows  containing  strong  shock  waves. 
To  assure  the  correct  shock  jump  relations  one  ought  either  to 
introduce  a  shock  fitting  scheme  or  else  to  use  a  difference 
scheme  in  conservation  form.   A  conservative  formulation  of  the 
small  disturbance  equation  has  been  given  by  Murman  [8],  and  the 
exact  potential  flow  equation  has  been  solved  in  conservation  form 
by  Jameson  [9]  for  flows  past  airfoils.   Comparisons  with 
experimental  data  show  no  clear  cut  advantage  to  using  the 


conservation  form  without  a  detailed  modeling   of  the  shock  wave 
boundary  layer  interaction  [10].   This  is  apparently  because  the 
error  in  the  shock  jump  relations  which  results  from  the  use  of 
the  nonconservative  schemes  is  in  the  same  sense  as  the  effect 
of  the  boundary  layer  interaction.  A  three  dimensional  scheme 
in  conservation  form  will  be  discussed  in  a  later  report. 


ANALYSIS 

Geometry 

Accurate  representation  of  the  finite  difference  boundary 
conditions  is  much  simplified  if  the  boundary  surfaces  lie  in 
coordinate  planes.   This  is  achieved  in  the  present  analysis  by 
a  sequence  of  transformations  based  upon  a  nearly  conformal  mapping 
of  the  physical  space  in  planes  containing  the  wing  sections, 
taken  in  the  streamwise  direction.   We  begin  by  considering  the 
physical  space  to  be  described  in  a  Cartesian  coordinate  system 
for  which  x,  y,  and  z  represent  the  streamwise,  vertical,  and 
spanwise  directions,  as  shown  in  Figure  1.   We  then  introduce  an 
arbitrary  singular  line,  just  inside  the  leading  edge  of  the 
profile  at  each  spanwise  station.   This  singular  line  will  be  the 
locus  of  branch  points  in  sxobsequent  transformations  in  each  of 
the  spanwise  planes  to  unwrap  the  wing  surface  to  a  shallow  bump; 
its  location  will  be  chosen  to  make  the  bump  as  smooth  as  possible. 
Representing  the  singular  line  as 

X  =  Xg(z) 


we  define 


y  =  ¥3(2) 


X  =  X  -  Xg  (z)  , 


y  =  y  -  ys(z)  /  (D 

z  =  z  . 

This  transformation  shears  out  the  wing  sweep  and  dihedral,  and 

puts  the  singular  line  at  the  origin  of  each  x,y  plane.   In  each 


of  these  planes  we  introduce  the  conformal  mapping 

(X,  +  iY.)^  =  2(x  +  iy)  , 


(2) 


which  maps  the  entire  wing  surface  to  a  shallow  bump  near  the 
plane  Y   =  0.   If  we  define  the  height  of  this  bump  as 


Y^  =  S(X,z)  , 


then  the  final  shearing  transformation 


X  =  X 


1  ' 


Y  =  Y^  -  S(X,z)  , 
Z  =  z   , 


(3) 


reduces  the  wing  surface  to  a  portion  of  the  plane  Y  =  0  • 

To  render  the  computational  domain  finite,  stretching 
transformations  are  introduced.   For  example, 

bY 


Y  = 


(1  -  Y^)^ 


0  <  a  <  1  , 


(4) 


is  used  to  map  the  planes  Y  =  +  <»toY  =  +  l.   Similar  transfor- 
mations are  used  outboard  of  the  wing  tip  in  the  Z  direction, 
and  downstream  of  the  trailing  edge  in  the  X  direction.   A 
sketch  of  the  resulting  rectangular  computational  domain  is 
shown  in  Figure  2. 

To  avoid  discontinuities  at  the  wjng  trailing  edge,  the 
branch  cut  in  each  spanwise  plane  is  continued  smoothly  down- 
stream.  In  the  physical  plane,  the  continuation  is  represented 


by 


In 


_    _* 
X  -  X 


_* 


y  =  Yte  +  ^(^te-^  ) 


_    _* 

X.  -  X 

te 


_   _* 

X  -  X 


(5) 


X.  -  X 

te 


where  t  is  the  mean  of  the  upper  and  lower  surface  slopes  at  the 

trailing  edge,   ^c.  ,  y    are  the  trailing  edge  coordinates,  and 

_* 

X   is  a  suitably  chosen  scaling  constant  (usually  taken  as  the 

ordinate  of  the  local  quarter-chord  point) .   In  the  solution, 
this  cut  is  taken  as  the    location  of  the  vortex  sheet,  across 
which  special  difference  formulas  must  be  applied.   Thus  we  make 
the  approximation  that  the  vortex  sheet  lies  in  a  fixed  surface 
necu:  the  plane  of  the  wing  which  leaves  the  trailing  edge 
smoothly  according  to  the  above  formula. 


Equation  of  Motion 

In  the  absence  of  strong  shock  waves,  the  steady,  inviscid 
motion  of  a  compressible  fluid  is  well  approximated  by  the  well 
known  equation  for  the  velocity  potential  <I>: 


(a^-u^)$  +  (a^-v^)$   +  (a^-w^)$   -  2uv«J)   -  2uw$   -  2vw<I>   =0,  (6) 
*       XX         yy  zz      xy      xz      yz    ' 


where  u,  v,  and  w  are  the  velocity  components  (i.e.,  the 
derivatives  of  $)  in  the  x,  y,  and  z  directions,  and  a   is  the 
speed  of  sound.   For  the  steady,  potential  flow  of  a  perfect  gas 
with  specfic  heat  ratio  y, 

2     2    Y-l  /  2^   2^   2,  ,_. 

a  =  a-  -  -'-2—  (u  +  V  +  w  )  ,  (7) 

where  a^  is  the  stagnation  speed  of  sound.   If  the  flow  is  uniform 
at  infinity,  parallel  to  the  x-y  plane,  and  inclined  at  an  angle  a 


to  the  X-axis,  the  far  field  singularity  can  be  removed  by  defin- 
ing the  reduced  potential  G   as 


G  =  $  -  X  cos  a  -  y  sin  a 

=  <t 


-  {7(^i-Yi)+Xg(z)|  cos  a  -  {xj^Y^+  y5(z)}  sin  a, 


(8) 


The  transformations  of  equations  (1),  (2),  and  (3)  applied  to 
equation  (6)  then  result  in  an  equation  of  the  form 


^  Sx  -^  «  Sy  ^  ^  Sz  -^  ^  Sy  -^  ^  Sz  -^  ^  ^yz  -^  ^  =  0  • 


(9) 


If  we  introduce  the  notation 


U  = 


V  = 


h  X. 


1  $ 
h  Y, 


^  =  -  ^1-  ^;  -  ^1-  y's ' 

X         y 

n  =    X,   X   -  Xi   y   , 
1-   s     i-  -^s  ' 

y      X 

r-  jx.cos  a  +  Y^s.in  a  +  G„-  S„G,,  }• 
h  (^  1  1  X    X  Yj 

7-  •l-Y^cos  a+  X,sin  a  +  G^  [ 


(10) 


(11) 


w  =  $, 


h^U  +  hnV  +  X  cos  a+  y  sin  a+  G„-  S_G^  , 

S  S  u  Ci      X. 


and 


where 


U  =  U  +  h^w  , 
V  =  V  +  hnw  , 


h2  = 


d(x  +  iy) 


d(X^  +  iY^) 


2     2 
=  ^1  ^  ^1  ' 


(12) 


(13) 


then  the  coefficients  in  equation  (9)  can  be  written  as 
A  =  a2{l  +  h^C^}  -  u2 

+  h^(a^-w^)S2  -  |2h^a^Cn  -  2Uv|s^ 

+  J2h^?a^-2hwulSj^S2  -  |2h^na^-  2hvjv\s^    , 

u2/  2     2\ 
C  =  h'|a   -w)-, 

D  =  -  2|a^(l+h^C^)-U^}Sj^+  |2h^Cna^-  2Uvl-  |2h^^a^-  2hwuls 


Z 


E  =  2h^Ca^  -  2hwU  , 


2   2   2 
F  =  -2h  (a  -w  )S„ 

Li 


(14) 


-  J2h^Ca^-  2hwu|sj^+  2h^na^-  2hwV  , 

R  =  l-ja^-d+h^?^)-  U^lsjjj^-  h^(a^-w^)S22-  |2h^U^-2hwulsj^2}^Y 

+  h3(a2-w2){{(x;2-y;2)x^__^  2x;y;x^__-  x^X^,-  y^X^_}u 
*•  ^  XX         xy      X      y* 

+  (-(x'2-y;2)x    +  2x;y;x,__+  x;x,_-  y;x^_}v} 
*•  xy         XX      y      x'  •' 

+  2h'^w{(X3^_x^-X^_y^)X^__+  (X^_x'+  X^_y^)  X^__|  (U^+V^) 
'■   X     y     XX     y      x     xy* 

+  i-  |x^U+Y^v|(U^+V^)  +  cos  a|h^(5^-n^)a^-U^-V^+h^(a^-w^)x^l 
+  sin  aJ2h^5na^-  2UV  +  h^(a^-w^)y^|  . 

Note  that  for  the  transformation  defined  by  equation  (2) , 

X,   =  X  /h^  , 

■^x    ^   ^  (15) 

Y^_  =  Y^/h^  , 

y 


and 

^12      2 

XX     h 

(16) 

Y 
12      2 
X,__  =    -i  (h^  -  4X^)  . 

The  siimnetry  condition  that  w  =  0  on  the  plane  z  =  0 
requires 

and  the  boundary  condition  that  the  flow  be  tangent  to  the  wing 
surface  requires 


y2 


^      (1  +  s^)  +  \s^  +   5S^  -  n|  K 


2  '      X'     \     Z  ^  X 


{-  -72  V  ^-^z-  ^^x-^^IK  ^  {-  ^z-  ^^x  ^  4^z 

j-X,_  cos  a-  X,_sin  afS  -  X.  cos  a+X-_sin  a  =  0  , 
^    X  V      -'       V        X 


(18) 


on  y  =  0, 

Downstream  of  a  finite  lifting  wing  there  will  be  a  vortex 
sheet.   Across  the  sheet  the  pressure  is  continuous,  but  there 
may  be  discontinuities  in  the  tangential  velocity  components. 
Convection  and  roll-up  of  the  vortex  sheet  are  ignored.   In 
reality,  the  component  of  velocity  normal  to  the  sheet  must   be 
zero,  but  in  our  approximation  it  is  simply  required  to  be  contin- 
uous.  Thus,  the  equation 

is  used  at  points  lying  on  the  vortex  sheet.  Also  the  disconti- 


nuity  in  potential  is  assumed  to  be  constant  along  streamwise 
coordinate  lines  downstream  of  the  trailing  edge.   The  value  of 
this  discontinuity  is  determined  by  the  Kutta  condition,  and  its 
spanwise  variation  determines  the  strength  of  the  vortex  sheet. 

Finite  Difference  Approximation 

The  success  of  the  type  dependent  difference  scheme  applied 
to  the  transonic  small  disturbance  equation  by  Murman  and  Cole  [1] 
can  be  attributed  to  the  fact  that  it  effectively  adds  a  direc- 
tional bias  to  the  equation  at  points  where  the  local  flow  is 
supersonic.  In  constructing  an  analogous  scheme  for  the  full 
potential  equation  in  general  curvilinear  coordinates  (which  may 
not  be  aligned,  even  approximately,  with  the  local  flow  direction), 
care  must  be  taken  to  ensure  that  this  bias  is  added  in  the  upwind 
direction,  i.e.,  in  the  direction  parallel  to  the  velocity  vector. 

A  method  with  this  property  has  been  proposed  by  Jameson  [6] 
To  illustrate  it,  we  return  to  the  potential  equation  in  the 
physical  coordinates.   The  equation  is  rearranged  as  if  it  were 
expressed  in  a  Cartesian  coordinate  system  aligned  with  the  local 
flow  direction,  s,  at  the  point  under  consideration.  Then 
equation  (6)  assiimes  the  canonical  form 

(a^-q^)$  ^  +  a^(V^$-$   )  =  0  (19) 

ss  ss 

where  q   is  the  magnitude  of  the  velocity. 


10 


The  relaxation  scheme  is  designed  to  simulate  an  artificial 
time  dependent  process  which  converges  to  the  desired  solution 
of  the  steady  state  equation.   In  the  finite  difference  approxima- 
tion to  the  potential  equation,  central  differences  are  used  to 
calculate  all  first  derivatives,  from  which  the  velocities  can  be 
determined  using  equations  (11) .   At  grid  points  where  the  flow 
is  subsonic,  central  differences  are  also  used  to  approximate  the 
second-order  derivatives  in  equation   (9) .  A  typical  central 
difference  formula  for  G  ^  is 

(n+1)       _    (n+1)  (n)       (n) 

G.  ,  .,-{-)  G.  .  ,  -  2(1-  -)  G.  .  ^  +  G._^,  .  , 
r        -      1-1, 3>k    CO    i,3>^<^ ^        i/]>k    i+l,3,k 

where  the  superscripts  denote  the  iteration  level   and  u  is 
the  relaxation  factor  [6].   If  we  regard  each  iteration  as  repre- 
senting an  advance  At  in  an  artificial  time  coordinate,  this 
formula  can  be  interpreted  as  an  approximation  to 

Sx  -  i  'St  *  s  <§  -  l'=t' 

Similarly,  the  formula 

p(n)         ^(n)         p(n+l)       ^  (n+1) 
_   _  ''i+l,j-H,k  "  ^i-H,j-l,k  "  ^i-l,j  +  l,k   "i-l.j-3,k     ,,,. 
^XY 4AXAY  ^"^^^ 

can  be  interpreted  as  an  approximation  to 

r   -  1  ^  r 
^XY    2  Ax  ^Yt  • 

The  relaxation  process  can  thus  be  regarded  as  an  approximation  to 
the  time  dependent  equation 


11 


(M^-l)G^^-  G^^-  G^^+  2a.G^.+  2a,G  .+  2a^G.+    6G.  =   Q  ,        (22) 
ss   nun   nn    1  st    2  mt    3  nt    t 


where  M  =  q/a  is  the  local  Mach  number,  m  and  n  are  suitably 

scaled  coordinates  in  the  plane  normal  to  the  velocity  vector, 

and  Q  contains  all  the  terms  in  the  equation  other  than  the 

principal  part.   The  coefficients  a, ,a_,a^,  and  6,  depend  on 

the  mix  of  old  and  updated  values  in  the  difference  equations 

as  well  as  any  explicit  time-like  or  mixed  terms  that  have  been 

added  for  stability. 

Introducing  the  new  time  coordinate 

""I 
T  =  t  -  —J —  s  +  a_m  +  a-n  , 

M'^-l       ^      -^ 
transforms  equation  (22)  to 

2 


(M^ 


In  order  to  ensure  the  convergence  of  the  scheme,  we  require 
that  equation  (23)  should  be  a  damped  three-dimensional  wave 
equation.   This  will  be  the  case  if 

a^  >  (M^-1)  {olI+olI)     .  (24) 

At  points  where  the  velocity  is  supersonic,  upwind  differ- 
ences are  used  to  represent  contributions  to  G   in  the  first 

ss 

term  of  equation  (19) .   This  is  done  using  formulas  of  the  type 
2Gf"-'\)  -  Gf").  ,  -  2Gf"rl).  ,  +  g!") 


r        -        i/J/k    i/j/k     i-l,j,k    i-2,j,k 

^(n+1)    ^(n+1)      ^(n+1)    .  ^(n+1) 
^   i,D,k    1-1, 3, k    i,:]-l,k    i-l,]-l,k 
XY  AX AY 

12 


(25) 


These  formulas  also  have  the  property  of  guaranteeing  diagonal 
dominance  for  the  updated  values  on  each  line.  The  formula  for 
G   can  be  interpreted  as  representing 

^XX  "^  ^  AX  ^Xt  • 
Together  with  analogous  formulas  for  G    and  G    ,  this  intro- 


duces  a  term,  equal  to 


2{M^-1)G  ^ 

St 


into  equation  (22) .   To  ensure  that  equation  (24)  is  satisfied 

2 
at  points  near  the  sonic  line  where  (M  -1)  is  small,  the  coeffi- 
cient of  G    can  be  further  augmented  by  adding  a  term  of  the 
form 


{"^Xt  ^  ^St  ^  h'^St}  ' 


6  ^  ^UG^^  +  VG„^  +  h"wG,^^  ,  (26) 


where   B  >  0  is  appropriately  chosen.   The  required    mixed 
derivatives  can  be  constructed  in  the  form 

g(n+l)  _  ^(n)    _  G^"'*'^)    +  G^^^ 
At     _   i,j,k     i/j,k    i-l,j,k    i-l,j,k  .^^. 

AX  ^Xt  ^P  •         ^      ' 

The  supersonic  difference  scheme  is  completed  by  using  central 
difference   formulas  similar  to  equations  (20)  and  (21)  to 
evaluate  contributions  to  the  second  term  of  equation  (19) , 
but  with  03   set  to  unity,  as  suggested  by  a  local  von  Neumann 
test  [6] . 
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Boundary  Conditions 

The  boundary  condition  at  infinity  is  particularly  simple 
because  the  square  root  transformation  reduces  the  entire  vortex 
wake  to  the  X-Z  plane  at  downstream  infinity.   Therefore,  since 
the  uniform  stream  singularity  has  been  removed  by  the  introduc- 
tion of  the  reduced  potential,  the  Dirichlet  condition 

G  =  0 
is  appropriate. 

On  the  X-Y  and  X-Z  planes,  finite  difference  approximations 
to  the  Neumann  boundary  conditions  specified  by  equations  (17)  and 
(18)  must  be  applied  to  those  portions  representing  solid  bounda- 
ries  (i.e.,  the  symmetry  plane  and  the  wing  surface).   At  the  wing 
surface,  central   difference  approximations  are  used  in  equation 
(18)  to  define  values  of  the  reduced  potential  at  image  points 
located  one  mesh  spacing  below  the  X-Z  plane.   A  similar  method 
is  used  on  the  symmetry  plane,  but  due  to  the  high  degree  of 
nonorthogonality   of  the  coordinate  system  when  the  wing  is  highly 
swept,  simple  central  differences  become  unstable.   Thus,  to  set 
the  potential  values  at  the  image  points  for  the  symmetry  plane, 
the  X-differences  required  in  equation  (17)  are  evaluated  by 
averaging  one-sided  differences  on  either  side  of  the  symmetry 
plane,  taken  in  the  upwind  direction  in  the  image  plane,  and  in 
the  downwind  direction  in  the  first  plane  in  the  flow  region. 
The  symmetry  condition  thus  remains  formally  second  order  accurate, 
and  the  incorporation  of  the  image  point  whose  value  is  being  set 
into  the  X-difference  adds  to  the  stability  of  the  scheme.  This 
method  of  handling  the  symmetry  condition  has  proved  statle  for 
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sweep  angles  in  excess  of   35  degrees. 

At  points  on  the  X-Z  plane  which  do  not  lie  on  the  wing 
surface,  the  values  of  the  reduced  potential  at  the  image  points 
are  taken  to  be  those  of  the  associated  point  on  the  other  side 
of  the  branch  cut,  allowing  for  a  discontinuity  across  the  vortex 
sheet.   The  value  of  this  discontinuity  is  taken  to  be  independent 
of  X  at  each  spanwise  station,  and  its  value  is  determined  by  the 
Kutta  condition  that  the  flow  leave   the  trailing  edge  smoothly. 

One  final  note  concerns  points   which  lie  on  the  contin- 
uation of  the  singular  line  outboard  of  the  wing  tip.   At 
these   points   the  mapping   is  singular,   and  a   special 
limiting  form  of  the  difference  equations  must  be  used.   At 
points  where  the  solution  is  regular,  the  nonlinear  terms  of 
the  potential  equation  are  of  0(l/h),  while  the  Laplacian 
transforms  to 


h^     Vi      Vi        2^ 

Thus,  in  the  limit  as  h  tends  to  zero. 


'11      11 
is  a  suitable  limiting  form. 


*X,X,  -^  *Y,Y,  =  °  (28) 
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RESULTS 
Computational  procedure 

The  potential  formulation  is  particularly  attractive  for 
three-dimensional  calculations  because  it  requires  the  storage  of 
only  one  quantity  at  each  grid  point,  and  the  number  of  grid 
points  required  to  accurately  describe  these  flow  fields  is  large. 
Even  so,  it  is  impractical  to  store  the  entire  solution  array 
in  the  high  speed  core  of  many  current  computing   machines . 
Fortunately,  since  the  analysis  presented  here  depends  on  a 
relaxation  solution  of  the  difference  equations,  it  is  not  neces- 
sary to  have  the  entire  solution  immediately  available  at  all 
times.   It  is,  therefore,  stored  on  a  disk  file,  and  read  into 
core  one  X-Y  plane  at  a  time.   At  any  time  during  the  solution 
procedure,  the  values  of  the  potential  on  four  such  planes  are 
in  the  core.   Old  values  are  buffered  in  and  new  values  buffered 
out  of  core  while  other  calculations  are  being  performed  as  much 
as  possible,  to  keep  the  process  efficient. 

In  each  X-Y  plane,  the  equations  are  solved  by  successive 
line  overrelaxation.   The  plane  is  divided  into  three  regions, 
as  shown  in  Figure  3.   In  the  central  region  the  equations  are 
relaxed  along  horizontal  lines,  sweeping  from  infinity  to  the 
wing  surface.   In  the  outer  regions  the  equations  are  relaxed 
along  vertical  lines,  sweeping  away  from  the  central  region  to 
infinity.   Such  a  sweep  pattern  ensures  that  the  sweep  direction 
will  not  be  opposed  to  the  flow  direction  in  any  supersonic  zones. 
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which  would  result  in  instability.   In  many  cases,  the  central 
region  can  be  taken  to  cover  the  entire  plane;  that  is,  only 
horizontal  line  relaxation  is  used. 

To  speed  convergence,  an  initial  calculation  is  usually 
performed  on  a  coarse  grid,  typically  containing  48x6x8  grid  cells 
in  the  X,  Y,  and  Z  directions  respectively.   This  solution  is 
then  interpolated  onto  a  finer  grid  containing  twice  as  many  mesh 
cells  in  each  direction,  and  is  used  as  a  starting  guess  for  an 
intermediate  solution.   The  process  is  repeated  once  again  to 
give  the  final  solution  on  a  grid  containing  192x24x32  mesh  cells. 
A  typical  run  consists  of  100  relaxation  sweeps  on  each  grid, 
requiring  a  total  of  approximately  85  minutes  of  CPU  time  on  a 
CDC  6600.   The  same  program  has  been  run  on  the  CDC  7600,  for 
which  a  similar  calculation  requires  about  15  minutes. 


Examples 

In  this  section  we  present  the  results  of  calculations 
using  the  swept  wing  program,  and  compare  the  predicted  surface 
pressure  distributions  with  those  measured  in  experiments.  The 
comparisons  are  made  for  two  different  wings,  each  typical  of 
a  class  of  swept  wings  of  the  subsonic  transport  type. 

The  first  wing  geometry  is  representative  of  the  tip 
panel  of  a  relatively  simple  wing  of  conventional  high  speed  sec- 
tion shape.  It  has  a  uniform  section  of  9.8  percent  thickness  ratio. 
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and  the  planform  has  a  leading  edge  sweep  angle  of  30",  a  taper 
ratio  of  0.7,  and  an  aspect  ratio  of  3.8.   A  program  generated 
projection  drawing  of  the  wing  is  shown  in  Figure  4.   The  wing 
was  tested  by  Monnerie  and  Charpin  [11]  of  the  ONERA,  and  carries 
their  designation  of  wing  M-6. 

The  first  results   presented  are  at  a  free  stream  Mach 
number  of  0.9226  and  zero  angle  of  attack,  resulting  in  zero  lift 
for  this  symmetrical  wing.  Figure  5  compares  the  calculated  and 
measured  streamwise   surface  pressure  distributions  at  the  20,  45, 
65,  and  95  percent  semispan   locations  [11,12],   Agreement  is 
quite  good,  including  the  predicted  shock  location. 

Figure  6  shows  similar  results  for  the  same  wing  at  a 
Mach  number  of  0.919  and  an  angle  of  attack  of  3.07  degrees. 
Again,  agreement  between  the  computed  and  experimental  results 
is  quite  good,  with  the  exception  of  the  shock  location  on  the 
lower  surface,  which  is  somewhat  further  aft  than  predicted  by 
the  calculation. 

Figure  7  shows  a  program  generated,  three-dimensional, 
projection  view  of  the  wing  surface  pressure  distribution  at  a 
Mach  number  of  0.840  and  an  angle  of  attack  of  3.06  degrees. 
This  is  a  particularly  interesting  case  because  of  the  merging 
of  two  shocks  into  one  on  the  wing  upper  surface  as  one  proceeds 
outboard.   This  pattern  is  graphically  illustrated  in  the  projec- 
tion view.   Figure  8  shows  comparisons  of  the  calculated  results 
with  experimental  data,  again  at  the  20,  45,  65,  and  95  percent 
semispan  stations.   Agreement  is  quite  good,  including  the 
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prediction  of  the  double-shock  pattern  at  the  inboard  stations. 

Figure  9  shows  the  projection  view  of  the  wing  surface 
pressure  distribution  at  a  Mach  number  of  0.837  and  an  angle 
of  attack  of  6.06  degrees.   Again,  the  calculation  predicts  the 
merging  of  a  double  shock  pattern  inboard  to  a  single  shock 
further  outboard.  Comparisons  with  data,  shown  in  Figure  10 
show  that  agreement  is  still  quite  good. 

The  second  geometry  is  representative  of  wings  being 
considered  for  the  next  generation  of  subsonic  transport  aircraft. 
The  wing  is  twisted,  both  aerodynamically  and  geometrically,  is 
highly  tapered,  and  has  a  discontinuity  in  trailing  edge  sweep 
angle  at  the  35  percent  semispan   location.   The  planform  has  a 
leading  edge  sweep  angle  of  35  degrees  and  an  aspect  ratio  of  7. 
It  has  5  degrees  of  dihedral.   It  is  defined  by  four  distinct 
streamwise  sections  (at  the  12,  35,  70,  and  100  percent  semispan 
stations) ,  with  linearly  interpolated  coordinates  between.  The 
streamwise  thickness  ratio  varies  from  16.3  percent  at  the  root 
to  11.9  percent  at  the  tip.   For  the  wind  tunnel  tests  the  wing 
was  mounted  on  a  quasicylindrical  fuselage  which  extended  to 
the  12  percent  semispan.   For  the  computations,  the  symmetry 
plane  was  assumed  to  be  at  the  same  spanwise  station  as  the 
wing-fuselage  intersection  in  the  tests.   A  projection  drawing 
of  the  wing  (extended  to  the  fuselage  centerline)   is  shown  in 
Figure  11.   For  these  calculations,  the  wing  geometry  was  modified 
to  account  for  boundary  layer  effects  by  adding  the  displacement 
thickness  obtained  from  two-dimensional  boundary  layer  calculations 
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multiplied  by  an  empirically  determ.ined  spanwise  weighting 
factor.   The  wing  was  one  of  several  tested  in  a  cooperative 
program  by  the  Douglas  Aircraft  Company  and  the  NASA  Ames 
Research  Center  in  the  Ames  11-foot  tunnel  at  a  Reynolds  number 
of  approximately  5x10  ,  based  on  the  mean  aerodynamic  chord. 

A  program  generated  three-dimensional  projection  drawing 
of  the  upper  and  lower  surface  pressure  distributions  for  this 
wing  is  shown  in  Figure  12.   (This  particular  case  was  run  with 
no  correction  for  boundary  layer  displacement  effect,  and  with 
the  wing  extended  to  the  fuselage  centerline.) 

Comparisons  with  experimental  data  are  shown  in  Figures 
13  and  14.   The  first  case.  Figure  13,  shows  streamwise  surface 
pressure  distributions  at  a  number  of  spanwise  stations  for  a 
Mach  number  of  0.75  and  an  angle  of  attack  of  2.2  degrees. 
Agreement  with  experiment  is  seen  to  be  excellent,  including 
the  location  and  strength  of  the  rather  strong  shock  near  the 
leading  edge  on  the  wing  upper  surface. 

Figure  14  shows  similar  comparisons  at  a  Mach  number  of 
0.84  and  an  angle  of  attack  of  1.8  5  degrees.   Again,  agreement 
is  quite  good,  although  the  resolution  of  the  first  (rather  weak) 
shock  of  the  inboard  double  shock  pattern  seems  lost  between  the 
35.5  and  50  percent  semispan  locations. 

The  results  displayed  in  Figures  13  and  14  were  kindly 
supplied  by  R.  M.  Hicks  and  P.  A.  Henne.   Further  details  of 
the  wing  geometry,  calculations,  and  test  conditions  are 
contained  in  [13] . 
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CONCLUSIONS 

A  numerical  method  has  been  presented  for  determining  the 
inviscid  transonic  flow  past  a  swept  wing.  The  method  is  based 
on  a  type-dependent,  finite  difference  approximation  to  the  full 
potential  equation,  solved  in  a  computational  domain  designed 
for  accurate  application  of  the  wing  surface  and  symmetry  plane 
boundary  conditions.   Calculated  surface  pressure  distributions 
agree  well  with  experimental  data  for  wings  of  conventional  and 
supercritical  section  shape  (when  the  geometry  in  the  latter 
cases  is  corrected  for  the  displacement  effect  of  the  boundary 
layer)  . 

Mapping  techniques  similar  to  those  used  here  could  be 
used  to  treat  more  realistic  geometries,  e.g.,  a  wing  mounted 
on  a  fuselage  [14].   The  recasting  of  the  finite  difference 
approximation  into  conservation  form  would  also  be  an  important 
theoretical  contribution. 

Finally,  as  was  mentioned  in  the  preceding  section,  these 
calculations  require  a  substantial  amount  of  computer  time. 
Thus,  methods  of  accelerating  the  convergence  of  the  iterative 
scheme  are  particularly  important  in  three-dimensional  problems. 
A  number  of  techniques  to  achieve  this  have  met  with  success  in 
two-dimensional  calculations,  including  a  hybrid  Poisson-solver/ 
relaxation  technique  [15,16],  a  multi-grid  method  [17],  and  an 
alternating-drection  method  [18] .   The  extension  of  these  methods 
to  three-dimensional  calculations  should  result  in  great  savings. 
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(a)   Plan  View 


(b)   Front  View 


Figure  1.   Geometry  of  Swept  Wing. 
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Figure   2.      Sketch  of   Computational  Domain. 
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Figure  3.   Sweep  Directions  in  Computational  Plane. 
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FIGURE  4.     GEOMETRY  OF  ONERA  WING. 
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FIGURE  7 
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FIGURE  9 


UPPER  SURFACE  PRESSURE 


LOWER  SURFACE  PRESSURE 
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FIGURE  11.     GEOMETRY  OF  DOUGLAS  WING. 
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FIGURE  12.     THREE-DIMENSIONAL  SURFACE  PRESSURE  DISTRIBUTION 
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Appendix  A.    Description  of  the  program 

All  the  numerical  results  in  this  report  were  generated 
by  the  computer  program  FLO  22  listed  in  Appendix  B.   This 
program  includes  options  to  treat  both  a  swept  wing  on  a 
wall  (Figure  Al) ,   and  an  isolated  yawed  wing  (Figure  A2) . 
For  swept  wing  calculations  the   sheared  parabolic  coordi- 
nates are  introduced  in  planes  parallel  to  the  free  stream. 
In  the  treatment  of  a  yawed  wing  the  whole  coordinate  system 
is  rotated  through  a  specified  yaw  angle,  so  that  the  X-Y 
planes  are  normal  to  the  leading  edge  of  the  wing  at  its 
center  line.   In  either  case  the  wing  section  can  be  varied 
in  an  arbitrary  manner,  and  the  only  restriction  on  the  plan- 
form  is  that  the  leading  edge  may  be  any  smooth  curve,  but 
it  should  not  have  kinks,  since  these  would  cause  the  second 
derivatives  of  the  singular  line  of  the  coordinate  system 
to  become  unbounded.   Kinks  are  permitted  in  the  trailing  edge, 
on  the  other  hand.   The  trailing  edge  defined  by  the  input 
is  actually  replaced  by  a  piecewise  straight  line  connecting 
the  nearest  mesh  points  in  the  computational  lattice. 
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The  geometry  is  defined  by  giving  the  wing  sections  at 
successive  span  stations  from  the  N'iug  root  to  the  tip,  or  in 
the  case  of  a  yawed  wing,  from  the  leading  to  the  trailing  tip. 
Up  to  11  span  stations  may  be  used  for  this  purpose,  and  the 
planform.   and  dihedral  are  determined  by  specifying  the  chord 
and  the  x  and  y  coordinates  of  the  leading  edge  at  these  span 
stations.   The  wing  section  at  each  station  is  then  determined 
by  scaling  and  rotating  a  prescribed  profile,  given  by  a  table 
of  X  and  y  coordinates.   If  the  wing  sections  are  similar,  only 
the  profile  for  the  first  station  need  be  read  in.   The  coordi- 
nates  for  the  other  stations  are  obtained  by  scaling  the  original 
profile  to  the  proper  chord,  and  rotating  it  to  obtain  the 
appropriate  twist.   If,  on  the  other  hand,  the  sections  are  not 
similar,  the  program  permits  the  coordinates  of  new  profiles  to 
be  read  in  at  each  span  station.  The  wing  section  between 
stations  is  generated  by  interpolation.  The  location  of  the 
singular  line  about  which  the  wing  is  unwrapped  by  the  square 
root  transformation  is  determined  by  the  parameters  XSING  and  YSING, 
which  must  be  specified  at  each  span  station.   It  is  important 
to  choose  these  so  that  the  mapped  profile  does  not  have  any  sharp 
bumps . 

The  main  input  to  the  program  is  read  from  Tape  5,  and 
the  output  is  written  on  Tape  6.   Tapes  1,  2  and  3  are  disk  files 
used  for  internal  storage  in  order  to  reduce  the  requirements  for 
high  speed  memory.   Tape  4  is  a  permanent  storage  device  such  as 
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a  magnetic  tape  en  which  an  intermediate  result  can  be  saved. 
The  computation  can  then  be  continued  for  more   iterations, 
starting  from  the  values  saved  on  Tape  4.   The  disk  instruc- 
tions in  the  version  of  the  code  listed  in  Appendix  B  are 
specialized  to  the  CDC  6600  using  the  FTN  compiler.   Otherwise 
the  code  should  be  readily  adaptable  to  other  computers. 
The  data  deck  for  a  run  is  arranged  to  include 
title  cards  listing  the  required  data  items.   The  complete  set 
of  title  cards  provides  a  list  of  all  the  data  which  must  be 
supplied,  and  can  be  used  as  a  guide  in  setting  up  a  data  deck. 
Each  title  card  is  followed  by  one  or  more  cards  supplying  the 
numerical  values  of  the  parameters  listed  on  the  title  card. 
All  data  items  are  read  as  floating  point  numbers  in  fields  of 
10  columns,  and  values  representing  integer  parameters  are 
converted  inside  the  program.   A  glossary  of  the  input  parameters 
is  given  in  Table  1,  and  a  typical  data  deck  is  shown  in  Table  2. 
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Table  1.  Glossary  of  input  parameters 

(Listed  in  order  of  their  occurrence  on  the  data  title  cards) 
TITLE  CARD  1 


NX 


NY 


NZ 
FPLOT 


XSCAL,  PSCAL 


FCONT 


The  number  of  mesh  cells  in  the  direction  of  the 
chord  used  at  the  start  of  the  calculation. 
NX  =  0  causes  termination  of  the  program. 

The  number  of  mesh  cells  in  the  direction  normal 
to  the  chord  and  span. 

The  number  of  mesh  cells  in  the  span  direction. 

Controls  generation  of  plots, 

FPLOT=0.  for  a  print  plot  but  no  Calcomp  plot 

at  each  span  station. 

FPL0T=1.  for  both  a  print  plot  and  a  Calcomp  plot 

at  each  span  station. 

FFL0T=2.  for  a  Calcomp  plot  but  no  print  plot  at 

each  span  station. 

FPL0T=3.  for  a  three  dimensional  Calcomp  plot  only. 

Control  the  scales  of  the  Calcomp  plots. 

XSCAL>0.  scales  each  section  plot  to  XSCAL 

XSCAL=0 .  scales  each  section  plot  to  5.0 

XSCAL<0.  scales  the  maximum  chord  to  XSCAL,  and 

each  section  plot  proportionately  to  the  local  chord, 

PSCALt^O.  sets  the  pressure  scale  to  PSCAL  per  inch 

in  each  section  plot. 

PSCAL=0 .  sets  the  pressure  scale  to  0.4  per  inch 

in  each  section  plot.    Also, 

PSCAL>^0.  scales  the  three  dimensional  plot  so 

that  the  span  or  semispan  is  5.   If  PSCAL=0.  and 

XSCALt^C.  then  the  three  dimensional  plot  is 

scaled  so  that   the  maximum  chord  is  1/2  XSCAL. 

Indicator  which  determines  the  manner  of  starting 
the  program. 

FCONT=0.  indicates  the  calculation  begins  at 
iteration  zero. 

FC0NT=1.  indicates  the  computation  is  to  be 
continued  from,  a  previous  calculation.   In  this 
case  the  values  of  the  velocity  potential  and  the 
circulation  are  read  from  a  magnetic  tape  where 
they  were  previously  stored  (Tape  4) .  It  is  still 
necessary  to  provide  the  complete  data  deck  to 
redefine  the  geometry.   The  count  of  the  iteration 
cycles  is  continued  from  the  final  count  of  the 
previous  calculation  and  the  maximum  number  of 
additional  iterations  to  be  performed  is  defined 
by  MIT. 
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TITLE  CARD  2 
MIT 

cov 


PI 


P2 


P3 


BETA 


STRIP 


FHALF 


The  maximum  number  of  iteration  cycles  which  will 
be  computed . 

The  desired  accuracy.   If  the  maximum  correction 
is  less  than  COV  the  calculation  terminates  or 
proceeds  to  a  finer  mesh,  otherwise  the  number 
of  cycles  set  by  MIT  are  completed. 

The  subsonic  relaxation  factor  for  the  velocity 
potential.   It   is  between  1.  and  2.  and  should 
be  increased  towards  2.  as  the  mesh  is  refined. 

The  supersonic  relaxation  factor  for  the  velocity 
potential.   It  is  not  greater  than  1.  and  is 
normally  set  to  1. 

The  relaxation  factor  for  the  circulation. 

It  is  usually  set  to  1.,  but  can  be  increased. 

The  damping  parameter  controlling  the  eimount  of 
added  (f  .  (see  equation  (2.6),  page  13). 

It  is  normally  set  between  0.  and  0.25. 

Determines  the  split  between  horizontal  and 
vertical  Line  relaxation  and  is  the  proportion 
of  the  total  mesh  in  which  horizontal  line  relaxa- 
tion is  used.   Fastest  convergence  is  usually 
obtained  by  setting  STRIP  =1.  so  that  horizontal 
line  relaxation  is  used  for  the  entire  mesh. 
If  convergence  difficulties  are  encountered  STRIP 
may  be  reduced  to  some  fraction  between  0 .  and  1 . 

Determines  whether  the  mesh  will  be  refined. 

FHALF=0.:  the  computation  terminates  after 

completing  the  prescribed  number  of  iteration  cycles 

or  after  convergence. 

FEKLF^O . :    the  mesh  spacing  will  be  halved  after  MIT 

cycles  have  been  run  on  the  crude  mesh  size.  An 

additional  data  card  must  be  provided  for  the 

refined  mesh  giving  the  numierical  values  requested 

by  Title  Card  2.  If 

FHALF<0   the  interpolated  potential  will  be 

smoothed  | FHALF |  times . 
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TITLE  CARD  3 

FMACH 

YAW 

ALPHA 


CDO 


TITLE  CARD  4 
ZSYM 


NO 


SWEEP 1 


SWEEP 2 


SWEEP 
DiHEDl 

DIEED2 


DIKED 


The  free  stream  Mach  number. 

The  yaw  angle  of  the  wing  in  degrees. 

The  angle  of  attack  in  degrees.    When  the  wing 
is  yawed,  ALPHA  is  measured  in  the  plane  normal 
to  the  leading  edge,  not  in  the  free  stream 
direction. 

The  estim.ated  parasite  drag  due  to  skin  friction 
and  separation.   It  is  added  to  the  pressure  drag 
(sum  of  vortex  drag  plus  wave  drag)  calculated 
by  the  prograiri  to  give  the  total  drag. 


Determ.ines  whether  to  treat  a  wing  on  a  wall  or 

an  isolated  wing. 

ZS5[:M=1.:  the  wing  is  on  a  wall 

ZSYM=0.:  the  wing  is  an  isolated  wing  at  a  yaw 

angle  given  by  YAW. 

The  number  of  span  stations  at  which  the  wing  section 

is  defined  on  subsequent  data  cards  from  the  wing 

root  to  the  tip   if  ZSYM=1.,  or  from  the  leading 

to  the  trailing  tip  if  ZSYM=0 .    If 

NC<3  it  is  assumed  that  the  wing  geometry  is 

the  same  as  for  the  last  case  calculated  and 

the  computation  for  new  values  of  FMACH,  YAW,  ALPHA 

and  CDO  begins  without  further  data  items 

being  read. 

Sweep  of  singular  line  at  the  wing  root  if  ZSYM=1., 
or  at  the  leading  tip  if  ZSYM=0 . 

Sweep  of  singular  line  at  the  tip. 
(SWEEPl  and  SWEEP2  are  used  as  end  conditions 
for  a  spline  fitting  the  x  coordinates  of  the 
singular  line.) 

Sweep  of  singular  line  in  the  far  field. 

Dihedral  of  singular  line  at  the  wing  root  if 
ZSYM=1.,  or  at  the  leading  tip   if  ZSYM=0 . 

Dihedral  of  singular  line  at  the  tip. 
(DIHEDl  and  DIHED2  are  used  as  end  conditions  for 
a  spline  fitting  the  y  coordinates  of  the  singular 
line. ) 

Dihedral  of  singular  line  in  the  far  field. 
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TITLE  CARD  5 
Z 

XLE,  YLE 
CHORD 

THICK 

ALPHA 


FSEC 


TITLE  CARD  6 
YSYM 


NU 
NL 

TITLE  CARD  7 
TRAIL 

SLOPT 


(The  geometry  at  the  first  span  station) 

Span  location  of  the  section. 

X  and  y  coordinates  of  the  leading  edge. 

The  local  chord  value  by  which  the  profile 
coordinates  are  scaled. 

Modifies  the  section  thickness.  The  y  coordi- 
nates are  multiplied  by  THICK. 

The  angle  through  which  the  section  is  rotated  to 
introduce  twist.  In  the  case  of  a  yav;ec  wing,  this 
angle  is  measured  in  the  axis  system  attached  to 
the  wing,  not  in  the  direction  of  the  free  stream. 

Indicates  whether  or  not  the  geometry  for  a  new 

profile  is  supplied. 

FSEC=0 . :  the  section  is  obtained  by  scaling 

the  profile  used  at  the  previous  span  section 

according  to  the  parameters  CHORD,  THICK,  ALPHA. 

No  further  cards  are  read  for  this  span  station, 

and  the  next  card  should  be  the  title  card  for  the 

next  span  station,  if  any. 

FSEC=1.:  the  coordinates  for  a  new  profile  are 

read  from  the  data  cards  which  follow. 

(Profile  Geometry  Supplied  if  FSEC=1.) 

Indicates  the  type  of  profile. 

YSYM=0.   denotes  a  cambered  profile.   Coordinates 

are  supplied  for  upper  and  lower  surfaces,  each 

ordered  from  nose  to  tail  with  the  leading  edge 

included  in  both  surfaces. 

YSYM=1.   denotes  a  symmetric  profile.  A  table 

of  coordinates  is  read  for  the  upper  surface  only. 

The  number  of  upper  surface  coordaintes. 

The  number  of  lower  surface  coordinates. 

For  YSYM=1.,  NL=NU  even  though  no  lower  surface 

coordinates  are  given. 

(Additional  Profile  Geom.etry  Supplied  if  FSEC=1.) 

The  included  angle  at  the  trailing  edge  in  degrees. 
The  profile  may  be  open,  in  which  case  it  is  the 
difference  in  angle  between  the  upper  and  lower 
surfaces. 

The  slope  of  the  mean  camber  line  at  the  trailing 
edge.   This  is  used  to  continue  the  coordinate 
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XSING,  YSING 


surface,  assumed  to  contain  the  vortex  sheet, 
smoothly  off  the  trailing  edge.   For   heavily  aft 
loaded  airfoils,  the  lift  is  sensitive  to  the 
value  of  this  parameter,  which  should  be  adjusted 
by  comparing  two  dimensional  calculations  using 
parabolic  coordinates  with  two  dimensional  calcula- 
tions in  the  circle  plane. 

The  coordaintes  of  the  singular  point  inside  the 
nose  about  which  the  square  root  transformation- 
is  applied  to  generate  parabolic  coordinates. 
This  point  should  be  located  as  symmetrically  as 
possible  between  the  upper  and  lower  surfaces  at 
a  distance  from  the  nose  roughly  proportional  to 
the  leading  edge  radius.   It  can  be  seen  whether 
the  location  has  been  correctly  chosen  by  inspect- 
ing the  coordinates  of  the  mapped  profile  printed 
in  the  output.   If  the  mapped  profile  has  a  bump 
at  the   center,  the  singular  point  should  be 
moved  closer  to  the  leading  edge.   If  the  mapped 
profile  is  not  symmetric  near  the  center,  with  a 
step  increase  in  y,  say,  as  x  increases  through  0, 
the  singular  point  should  be  moved  closer  to  the 
upper  surface.   The  coordinates  of  the  singular 
point  are  chosen  relative  to  the  profile  coordinates 
supplied  on  the  cards  which  follow. 


TITLE  CARD  8    (Upper  Surface  Coordinates) 

X,Y  The  coordinates  of  the  upper  surface.   These  are 

read  on  the  data  cards  which  follow,  one  pair  of 
coordinates  per  card  in  the  first  two  fields  of  10, 
from  leading  to  trailing  edge  inclusive. 

TITLE  CARD  9    (Lower  Surface  Coordinates,  Read  if  ISYM  =0.) 

X,Y  The  coordinates  of  the  lower  surface,  read  from 

leading  edge  to  trailing  edge.  The  leading  edge 
point  is  the  same  as  the  upper  surface  leading  edge 
point.   The  trailing  edge  point  may  be  different  if 
the  profile  has  an  open  tail. 

TITLE  CARD  10,11...  (Geometry  at  the  Other  Span  Stations) 

These  title  cards  are  the  same  as  Title  Card  5 
(geometry  for  the  first  span  station) .  The  number 
of  such  cards   depends  on  the  number  of  input  span 
stations  NC .   If  the  profiles  are  similar  at  each 
station  except  for  scaling,  thickness  to  chord  ratio 
and  rotation  to  introduce  twist,   FSEC=0.  and  no 
new  profile  coordinates   are  needed. 
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TABLE  2. 

DATA  DECK  FOR  ONERA  M6  WING 

""\^  Columns 
Cards  ^\.„^^^ 

1-10 

11-20 

21-30 

31-40  1  41-50 

1 
1 
I 

51-60 

61-70 

71-80  i 

Title  of  case 

ONERA  M6 

WING  (copied  onto 

output  aind  Calcom^  plots)  j 

1     .    i        ' 

1 

Title  Card 

NY 
48. 

NY 
6. 

NZ 
8. 

1        t 
FPLOT  '  XSCAL  [ 

1.      1    0.      : 

1         i 

1 
PSCAL   ' 

0.     ' 

FCONT 
0. 

1 

Title  Card 

MIT 
100. 
100. 
100. 

GOV 
l.E-6 
l.E-6 
l.E-6 

PI 
1.6 
1.6 
1.6 

P2       1  P3       ; 
i              1 
1.   ;  1. 

1.      i    1. 
1.      -1. 

BETA   i 
.10 
.10  . 
.10 

STRIP 
1. 
1. 
1. 

FHALF 
1. 
1. 

Title  Card 

MACH 
.840 

YAW 
0. 

ALPHA 
3.06 

CDO                     ,                ! 

i 

.010                   ' 

Title  Card 

ZSYM 
1. 

NC 
6. 

SWEEP 1 
29.9 

SWEEP 2    SWEEP 
29.9     29.9 

DIHEDl 
0. 

DIHED2 
0. 

DIHED 
0. 

Title  Card 

Z 
0. 

XLE 
0. 

YLE 
0. 

CHORD    THICK 
.6737    1. 

ALPHA 
0, 

FSEC 
1. 

Title  Card 

YSYM 
1. 

NU 
72. 

NL 
72. 

1 

1 
1 

1 

Title  Card 

TRAIL 
7.06 

SLOPT 
0 

XSING 
.00725 

YsiNG    :             ; 

0.        ■            i            !            !           1 
1                                      ' 

Title  Card 
(72  cards) 

X 
(Coordin 

Y 
ites  of  p 

rofile) 

j     (Upp 

1 
er  Surface)      j 

i 

[ 

Title  Card 

Z 
.2 

XLE 
.1150 

1  YLE 
0. 

CHORD 
.6147 

THICK 

1. 

ALPHA 
0. 

FSEC 
0. 

Title  Card 

Z 

.4 

XLE 
.2300 

YLE 
0. 

CHORD 
.5558 

THICK 
1. 

ALPHA 
0. 

FSEC 

0. 

i 

Title  Card 

Z 
.6 

XLE 
.3450 

YLE 
0. 

CHORD 
.4968 

THICK 
1. 

ALPHA 
0. 

FSEC 
0. 

Title  Card 

Z 
.8 

XLE 
.4600 

YLE 
0. 

CHORD 
.4379 

THICK 
1. 

ALPHA 
0. 

FSEC 
0. 

Title  Card 

Z 
1.0 

XLE 
.5750 

YLE 
0. 

CHORD 
.3789 

THICK 
1. 

ALPHA 

0. 

1 

FSEC 
0. 
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Both  graphical  and  printed  output  are  provided.   The 
wing  sections  defining  the  geometric  configurations  are 
printed  for  each  span  station,  if  they  are  different,  or  for  the 
first  span  station  only  if  the  sections  are  all  similar.   The 
program  next  prints  the  coordinates  of  the  unfolded  sections 
produced  by  the  square  root  transformations  at  the  root  and 
the  tip.   These  should  be  inspected  to  see  that  they  are  reason- 
ably smooth.   The  prograiri  alsc  prints  a  chart  of  an  indicator  IV 
showing  the  configuration  of  the  wing  in  the  coordinate  surface 
to  which  it  has  been  mapped.   The  values  of  IV  are  as  follows: 

IV  =  2   indicates  a  point  on  the  wing 

1   indicates  a  point  on  the  trailing  vortex  sheet 
0   indicates  a  point  on  the  singular  line 
-1   indicates  a  point  adjacent  to  the  edge  of  the  wing 

or  vortex  sheet 
-2   indicates  an  ordinary  point  not  in  contact  with  the 
wing  or  vortex  sheet. 

The  program  next  displays  the  iteration  history.  The 
maximum  correction   to  the  velocity  potential  and  the  maximiom 
residual  of  the  difference  equations  are  printed  at  each  cycle, 
together  with  the  locations  of  the  points  where  these  occur 
in  the  computational  lattice,  and  also  the  relaxation  factors, 
the  circulation  at  the  wing  center  line,  and  the  number  of 
supersonic  points. 
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After  a  specified  maximtun  number  of  cycles  has  been 
completed,  or  a  convergence  criterion  has  been  satisfied,  the 
section  lift,  drag  and  moment  coefficients  are  printed  for 
each  span  station,  and  the  pressure  distribution  is  printed 
or  displayed  in  a  Calcomp  plot  as  desired.   Finally  the  charac- 
teristics of  the  complete  wing  are  printed.  These  include  the 
coefficients  of  lift  and  form  drag  computed  by  integrating 
the  surface  pressure,  and  the  ratio  of  lift  to  form  drag. 
An  estimate  of  the  friction  drag  coefficient  may  be  supplied 
in  the  input,  and  this  v;ill  be  included  to  provide  an  estimate 
of  the  total  drag  coefficient  of  the  ratio  of  lift  to   total  drag, 
The  pitching,  rolling  and  yawing  moments  are  also  computed 
and  printed.   In  the  case  of  a  yawed  wing  these  are  in  an  axis 
system  normal  to  the  wing  leading  edge  at  its  center  line.  In 
the  case  of  a  wing  on  a  wall  the  rolling  mom.ent  is  the  root 
bending  moment. 

Finally  additional  Calcomp  plots  are  generated  if  they 
are  desired.  These  show  the  convergence  history,  and  also  a  viev; 
of  the  complete  wing  and  the  three  dimensional  pressure  distri- 
bution over  the  upper  and  lower  surfaces  separately,  with  the 
wing  root  or  the  leading  tip  at  the  bottom  of  the  picture.  If 
the  mesh  is  to  be  refined  the  program,  then  completes  the  sam.e 
sequence  of  calculations  and  output  for  the  new  mesh. 
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rection 
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Figure  Al.   Swept  wing  on  a  wall 
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Flow  direction 
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Figure  A2 .   Yawed  wing, 
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APPENDIX  B 


LISTING  OF  THE  PROGRAM 


PROGRAM  t^LOZ? 
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4 
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JC 
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^1 
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CJOKIN 

tSG^#^•A 
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),32(2b 
# v2(35) 
),xZZ(i 
IfKTt  Z, 
*>*  ALFHA 
P3,Bf TA 
S(2^1»l 
i.).  Yl  [  ( 
), l3( 11 
't  1  )  .  C  1  ( 
)»bV(19 
L{3i)  ,S 
3),F^C( 
aLF(3  ), 


i,TAPt2,  TAPt  3*TAPE't, 
"UUTPlT) 

IN    TKAiNSQMC    FLOW 
ATES 

kCh    197<i 

ASTCNY  JA1ESQN»JEC  197i-DtC  1976 

J5),cO(131),ZO(131) , 

f  iTcdOS), 

193), A3{193), 

) .B3(26), 

»03(  3:>), 

S), YC( 3t ),YZ(35) , YZZ  (35), 

iS>'f  ,KSY1,SCAL,SCALZ, 

,CA,5A,FlACH,M,N2,.>l3,ia 

, FR,IR,JR,KR,0G,1G, J&,KG,NS 

i), 

il),SLQPT (Ii),TPAIL(ll),NP(ii), 

)  ,  F  ^  (  1 1  )  , !:  5  (  1 1 )  , 

2  <.i),D2(24l), 03(2^.1), 

3),Sf1(  193),CP(i93), 

C0(35),SC'1(35),TITLE(20), 

3),P2C(3),P30(3)»aETA0(3), 

kES(  5o1),CUiJNT(501) 


7,c 

» cC 
,2) 
GGP 
REE 
ING 
5  3C 
630 
50t 
^IC 
NX 


9'j779i>13wd23 
0) 


hanTQNY  JAMtSON,COURANT  INSTITUTE/ 
NG  ANALYSIS  IN  TRANSCnIC  FLOW, 
Lie  CGiJ(«DlNATtS) 


AM  FLQa2,7CX,32 

DIMEnSIGNAL  WI 

SHEaRcU  PAkABC 
)  TITLE 
)  TITLE 
) 
)   R,x,FNY,FNZ,FPLQT,XSLAL,PSCAL,FCCNT/Fa1 


51 


NY         «  FNY 
NZ         •  FNZ 
If-  (NX.LT.l)  GD  Id  3Cx 
KPLQT      «  ABS(FPLCT) 
RtAO   (IkEAO»t>OC) 
Nh  «  0 

11  NM  »  NM   ♦! 

PtAO       (IRI:A0,5K)     F  I  T  (  SI )»  CO  VO  (  N"  ),  Pi  ^  (  N'i ),  ?<iO  (  N  ^1 )»  P3  J  (  Nn  )  * 
1  -^lTAOCNM),  S7I-  IPO(N^),  FHALF(Nf^) 

IF     (FHALF(Nil)  .Nt.C.  ,/sNO.Nrt.Ll  .2  )     GU    TJ    11 
FHALFO)       «    C. 
RtAD       (IPtA0,5OC) 

kEAO       (IRbA0»51C)     t-MACHf  rA>  AL,Ci)C 
YAW  =YA/kAD 

ALPHA  «     aL/PAO 

CALL    GECM       ( MD, NC» NP, ZS » XS* YS> X ll , YL t » 5 L DF 1 » TkA IL, X P, YP, 

1  SWFtPl,SWi.cPt»5»^H.  P,0lHE01»LIHeC2>OHhLi, 

2  XltC*CHGKD>^,ZTlP*  iSY10#Ki)Yr) 

ISYM      »  iSYro 

IF  (^LPHA.NE.C. )  l>Yh  ■  0 
IF  (KSYf  .Nt.O)  YAK  «=  0. 


CYAW 
SYAw 
CA 
SA 


CLS( YAa) 
SIN( YAW) 
CYAw*Cub (ALPHA) 
C YAU*SIN( ALPHA) 


62 


91 


IF  (FCONT.LT.l. )  GG  l L  vl 

READ   (<»)    NX*NY/NZ,N^/Kl,K2#Mr 

.MX  «  NX   +1 

MY  ■  NY   +E 

hZ  »  NZ   +3 

DC  62  K«1,MZ 

READ   Ci)  (  (G(  l.Jf  i)#I'i*MX  ),  J  =  1,MY) 

RUFFlR  QUT(N3#1)  <G(lf 1»1 )f G(MX,MYf 1)  ) 

IF  (LMT(N3).GT.G.)  Gu  TL  i 

bUFFtR  QUT{N1»1)  (  G  (  1 ,  1»  1  )  »  G  (  MX  ,  (^  Y,  1 )  ) 

IF  (UNIT(M)  .Gf  .0.  )  bU  TC  1 

CCNTINUE 

READ   Cf)    ( E0(K),k=k1,k2) 

PEwlND  N3 

fit.^iND  Nl 

f-EWlND  A 

CALL  COaKC 


(NX,NY»NZ»t<>bYr»XU(,ZTIP,XM^X,ZnAX, 
SY»iiCAL#:^CALZ»AX*AY»AZ» 
AI»A1.«2»«3*60*B1,62»B3,Z»^1.C2,C3) 
(NC#NZ#KSYM,K7ti*KTf2*CHaRUC» 
SwEEPl»SWEbP2»bwEEP»0IHEDl,0lHED2/01HtP» 
Zi,XLl:,YLL>xC»XZ,XZZ»YC*YZ>YZZ» 
Z>Ci*C2»C3#c:l»t2#F3*£^*E5*lN0) 
(ND,Nt»NC».'iX>NZ»lbY^i>KSYM»KTEl#KTt2»i>CAL» 
YAW*  AL»Z»Zj*XC»YC#bLa?T#  TRAIL/ XS*Yi)»NP. 
iTEl,]Tt2/iV»SO,ZC/AP,YP»01/02»D3/X,Y,  INO) 
IF  (IND.tU.O)  GC  Ui  2vi 
IF  (FCONT.GE.l. )  GU  TC  101 
NM  -1 


CALL  SINGL 


CALL  SURF 
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101 


103 

106 
11? 

ll<t 

lie 

116 

120 
122 


MT 

CALL 

IF  (I 

REwiN 

REWIN 

WkIT«^ 

FC3NT 

MIT 

KIT 

II-  (N 

dIT 

K»«ES 

JPEb 

CUV 

STRIP 

PlTA 

MX 

MY 

HZ 

KT 

Kl 

K2 

IF  (K 

Kl 

LZ 

If     (K 

WRITE 

FL.R^A 


f.GT. 


1 
2 


DC  10 
W^ITE 
WRITE 
WRITE 

FCRMA 


1 
2 


Df.'  11 
*RITt 
WRITE 
FCiPMA 
WkITt 
WRITE 
WRITE 
FCRiA 

DC  12 
WRITE 
WRITE 
FGR1A 
WRITE 
WRITE 


ESTIM 

Q.EC.O) 

D  N3 

D  Nl 
( IwRITf 
0 
F 

A 
N 
( 

C 

c 
c 
s 

e 

N 
N 
N 
N 
2 
N 

S  Y  M  .  c  0  .  > 

*  3 

■  S 

«  N 

SYM.Nt.D 

(  IwkIT, 

TC-^BHOIN 

27H  IN 

27hG( ( 

C  I«2»NX 

(IWHT, 

(IWKIT, 

(IWPIT* 

1 (49HCCH 

3^H  AN 

15H0 

^  I»2»NX 

(IWRIT* 

(IWPIT» 

TdbhO 

(IwRIT, 

(1*RIT, 

(IwRIT» 

T(H6H0Na 

15H0 
0  J«2»KY 
(I*RIT» 
(IWRIT, 
T(1&H0 
(IwRIT, 
(IfcRIT, 


GC  Tu  1 


tCi) 

t 

I T  ( N  M )  +  N  n 

IT 

^C.FHALF(NM),EU.C. )  KIT  ■  10 

IT 

"11   -MT  -dUt^OQ       *Z 


CVC(NM) 
TClPO(Nn) 
ETAO(NM) 
X   +1 

Y  ♦? 
Z  *i 

Y  +1 


)  GO  TO  iu3 

Z  *Z 

)  IZ  «  J 
It*) 

LICATK-N  uF  L'JCATILN  OF  wING  AND  VORTtX  SHEcT* 

CCGkUlNAft  PLANE  Y  «  0./ 
IV(I»K)jK«Kl,K^),I«2»NX)) 

bbi   )  (  IV(  I»r  )»K«K1,K2) 
tlC) 
112) 

ORCwIbb  CtLL  CISTRlbUTILN  IN  SQUARE  RCOT  PLANt* 
I  rAPPtu  ^.URFACL  CUCRDINATES  AT  CENTtR  LINE  AND  TIP/ 
X       ,i5H    RCGT  PP0FILE»15H    TIP  PROFILE  ) 

61C)  AC(I)»SOII*LZ)*SO( I»KTE2) 
lit) 

TE  LGCaTICn  »15H     PQwER  LAW   ) 
61C  )  XMAXfAX 
60C) 
lit) 

RMAL  CtLL  OlSTRIbUTIGN  IN  SQbARfc  ROOT  PLANE/ 
Y        ) 


eic)  BC(J) 

122) 

SCALF  FACT0R#i5H 
1 1 C  )  S  Y  *  A  Y 
60C  ) 


PGWcR  LAW   ) 
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12A 


WRITt    (IWRIT»12A) 
FORrlAT(^5hOSPANWiSt 


Z 

xZ 
YZZ 


1  lt>HO 

2  15h 

3  15H 
DO    126    K«K1»K2 

126    WkITt     (IWI<IT>61C) 

WRITE     (IWfIT,12fc) 
128    FGRiATd&HO       TIP    L0CATiQN*15H 

WRITE     (IWRIT,610)     ZMaX#AZ 

WRITt    (lWRlT#60t) 

WRITE    (IwPIT,132) 
132    FtjRiAT(19H0ITERATIVE    JULbTlCh/ 

1  ^jHOilflf    WIDTH    hLt<    hDRIZQMAL 

WRITE    (IWkIT,61C)    STkIP 

WRITf^    (IWPITflS't) 
i3<»    FuRflATdbHO 

WRITE     (IwHT,6^0) 

CALL    SECOND(T) 

WRITE     (IWRIT»70C)    T 

WRITE     (IWKlT,13t) 
136    FCRMATdiJHO  rACH    NO  .  IbH 

WRITE     (IwRlTfblC)     FMACH,TA,AL 

WRITE     (l6RIT»136) 
138    FOR'iAT(10HOITEFATiON,l3H  COkRtCTlGN 

1  I'J^  KLSIDUAl 

2  lOh  CIRCULATN,10H  RtL  FCT  IdCH 


CtLL  DISTRIBUTION  ANO  SINGULA-^  LiNt/ 

,15H       X  SING  »i5H       Y  SING 

>i:?H          YZ  »15H         Xll 
) 


Z{K),XC(K),YC(K),XZ(K),YZ(K)»XZZ(K),YZZ(K) 


PQfliK    LAW   ) 


Ll-^fe    RELAXATION) 


NX  »i: 

NX,NY»NZ 


NY 


15H 


Hi 


YAW 


,15H        AN&    OF     ATTACK) 


REL 


#4H 


FCT    2»xOh    REL 


K 
FCT 


i» 


1^1 


\  lUH  fit]  A 

NIT  =    NIT       +1 

JIT  •    JIT       +1 

PI  «=    PiO(Nf) 

P2  •    P2C(NM) 

P3  »    P3C(N»^) 

IF    (NlT.LE.iC)    PI    « 
IF     (MT.Lt.lO)    F3    » 
CALL    1IXFL0 
IF    (IL.PC.O) 


>iOH    SbMC    PIS) 


GC    TO    IM 


JO 

REWIND 

Rewind 

N 

Ni 

N2 

N3 

WRITt 


0 


Nl 

N2 


»    Nl 
«    N2 
»    N2 
■    N 
(iWHT,660) 


NlT#L;G*IG»j6»K&#FR>IPfJR*KR#tO{LZ), 
Plf P2»P3»BFTA>NS 
■H 

jRfcS  =  1 
TO  1^3 
♦  1 
-i 


1^3 


JRES       ■  JRES 
IF  ( JRES.EC.KRES) 
IF  (JRFS.NE.I)  GO 

NREs  ■  nrf:, 

CLiUNlT(NRtS)    =    Ml 

RtS(NReS)     «    FR 

IF     (JIT. EC. KIT)     GG    TO    251 

IF    (MT.LT.KIT.AND.Ati(DG) .GT.CQV.AND.A6S(0G)  .LT .10.)    GC 


TG    i<tl 
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l5l 


161 


60  TC  161 
IF  (jC.tO.l) 

ktWiND  Nl 

REWIND  N2 

JO 

N 

N3 

N2 

Nl 

GO  TO  1^1 

Rate 


GC  TO  1 


1 

N3 
N2 
Nl 

N 

0. 


IF  (NRES.GT.l)  RATt 


WRITE  (IwRIT,16i;) 
162  FCRMATdbHO  f^AX  KtS 
1         15H   kECLCTN 
Whirt  (H.K1T,67C)  k 

Call  second(T) 

WRITE  (IwRIT#70C)  T 
WRITE  (IWR1T,6CC) 
DO  16<r  L-l#3 
9LFFER  IN  (Nl,l)  (G 
IF  (UNIT(M).GT.O.) 
CONTINUE 


«  (  ABS(t<ES(Nkti) /KES(l)  )  ) 

♦♦(1./(CCUNTINRES)   -CJLNT(l))) 

IDAL  i»15H   f.AX  RbSIDAL  2,15H 

/CYCLE) 

Ei,(  x)»K:S(^PLS)»COU^T(^RES)»RATE 


WORK 


(l*i»L ),G(MX,rY*L) ) 
60    TO    lt>l 


16A 


171 


LX 

K 

K 

IF 

DC 

DC 


172 


«    NX/2      ♦ 
■    2 

'    K       ■*! 
(K.EU.NZ)  GO  TO 
172  J-1*MY 
172  I=1»MX 
C.  (I  >  J .  1  )   »  G  ( 1 ,  J  i  2 
Gd,  J»2)   «=  Gd,  J#3 
BUFFER  IN  (Nl,l)  (G 
IF  (UNIT(M)  .GT.O.) 
IF  (K.LT.KTEl.CR.K. 

11  «  ITEKK) 

12  «    ITE2(K) 
CALL    VELQ    (K,2»SV#S 
CHORD(K)        -    X(I1) 
CALL    FOPCF     (I1>U*X 
IF    (kPLQT.GT.  LAND. 
WRITE     (1WRIT,6CC) 
*RITC     (IWRIT»162) 
FCS*1AT(24H0StCTi0N 

1  15H0  fiACH 

*.ftITE    (IWRIT»61C)    F 
WRITE    (IWRIT»184) 

184  FQRMAT(15H0       bPAN    S 
1  15H  C 

185  WRITE  (IWRIT#610)  Z 
IF  (KPLOl.LE.l)  CAL 
IF  (KPLOT.LT.l.OR.K 
CALL    GPAPH    (IPLrT»I 

1  Z(K),SC 


l<il 


) 

(1*1#3)*G(MX,MY»3)  ) 

GO   TO   i;>i 
G1.KTE2)    GC    TC    171 


M>CP»X,  r) 

-X(LX) 

, Y*CP#ALf CHORD (K),XC( K),iCL(K ),  SCO(K),SCM(K)) 

K.GT.KTEl)  GC  TO  185 


182 


CHARACTtklSTICS/ 

NO    #15H         YA* 
HACri»YA» AL 


CL 


,15H   ANG  OF  ATTACK) 


*15ri 


TATiGN,l&H 

H       ) 

(K)»SCL(K)>SCD(K),SCM(K) 

L  CPLuT  (I1*I2*FMACH,X, Y,CP) 

PLaT.GT.2)  GO  TO  i71 

i,i2»X,Y,CP,TITLE»Fr'ACH,YA,AL, 

L(K)*  SCl>(K  )*ChORDO#XSCAL*PSCAl> 


CO 


55 


19] 


GC 
CAL 


TD  171 
L  TOTFCR 


CUl 

CD 

VLD 

!»- 

VLD 

IF 

WRI 

kRI 

192    FOR 
1 
WRI 
ViKI 

19<r  FGR 
1 
WRI 
WRI 
FOR 
WRI 

IF 

CAL 

CAL 


196 


(ABS(CDi 

{APS(Cl  ) 
Tf  (I*kl 
TE  (IWRI 
IAKZIHO 

IbtiQ 
TE  (I.<kl 
TE  (I*RI 
MAT(15mO 

15h 
TE  (IWPl 
TE  (I WRI 
MAT(15hO 
Tt  (I».'RI 
IND  Nl 
(KPLOT.L 
L  RPLCIK 
L    THREED 


{KTEi»KTti;»CHi:i^D»SCL*SC:#iCf»Z»XC» 

CL*Cu1»CMH*  CMK,CrY) 

CY4w*CDl 

CDO       +CDi 

0. 
).GT.i.L-t)     VLux    =    CL/COi 

0. 
.GT.l.E-6)     VLD    «    CL/CD 
T,60C) 
T, 192) 
WING    CHARACT£^ISTICS/ 

rACH    rj  ,li»H 

T»61C)     F^'ACrl,  YAjflL 
T»19'.) 

CL  filjH 

CD  »i5H 

T,6lL)    CL»CDi,CCO,CCi VLC1>VLD 
T,196) 

CM    PITCH       ,15H  C*1    ><CLL 

T»61t)    C!^PiCMR#CMY 


YAW 


CD    FCKM 
L/D    FaR^'. 


♦  lt?H       AhG    LF    ATTACK) 


>15H         CD    FkiCTlDN    » 
*15H  L/D  ) 


15H        CM  YAw 


(IC.EC.O 
(ISTOP.c 
(FHALF(N 


IF 
201  IF 
IF 
NX 
NY 
NZ 
CALL  COORD 


1 
2 

( 
1 
2 

3 

( 
1 
? 

IF 

CAL 

IF 

RtW 

KEW 

NSM 

IF 

DO 

CAL 

IF 

c\jZ    RfW 


CALL  SINGL 


CALL  SURF 


T.  1)  G 
IPLCl, 

(iPLcr 

VLD,C 

)  GC  T 

0.1)  G 

M) .EC. 

NX   ♦ 

NY   ♦ 

NZ   ♦ 

(NX,  NY 

SY,SC 

A0»A1 

(NC»NZ 

Sweep 

ZS#>L 

z*ci, 

(ND,Nt 
YAW,  A 
ITEl, 

0)  GC 


0  TO 
NRES* 

L  ♦  C  L'  * 
Q  151 
0  TO 
0.)  G 

NX 

NY 

NZ 

fHlfK 

AL,SC 

,A2,  A 

1,  SWt 
£,YLE 
C2»C3 

,NC,N 

ITE2, 
TO  c9 


( IND.tt. 

L  RtFIN 

(IG.E&.O)  GO  TG  'dZl 

iND  Nl 

INU  N2 

nr 

(NSMCG.L 

2C2  N-1, 

L  SMOO 

(IC.tU.C)  GO  7  3  221 

iND  Nl 

I K  Q  K'  2 


^01 

ktS*COUNT#TITL£,FMACH» YA, AL, NX,NY*NZ) 
H,Cf'*X#  Y*TITl£,  YA,AL* 
CHCRDC,XSCAL,PSCAL) 

0    TO    1 


SYM,XTEO*ZTIP»XMAX,Zf^AX, 
ALZ,AX#AY,AZ, 

i,bO,Bl,E2,b3,Z,Cl»C2,C3) 
*Kltl,HTt2f CHGkDG, 
EP2o»*EfcF,L;lHECl,DlHEDt,DlHED, 
,XC,XZ,XZZ,YC. YZ,YZZ» 
,tl,c2,f  :J,t^,Et;,  INO) 
A,NZ*ISYM,KSYi-1,KTtl#KTt2,SCAL, 
i. XCYCSL  OPT, TRAIL, XS»Y3#NP, 
iV,S0,ZL,»XP,rP,01,Di:,03,X,Y,  INC) 
i 


-FHAlFCN^) 
T.x)    on    lu    211 

NSMCD 


56 


211  N  ■  M 

M  "  N2 

N2  »  N3 

N3  ■  N 

NM  »  NM 

MT  «  0 

GO  rc  101 

221  NX 

NY 
N2 
CALL  COGRL 


♦  1 


=  NX/2 

•  NY/2 
«  NZ/2 
(NX,^rfNZf^br^'»XTfO»2TiP»XMAX»2MAX> 

1  SY,SCAL» itALZ» AX, AY,A2» 

2  AO»Al#A2#A3,6J,Bl,f-2*e3#Z,Cl»C2»C3) 
CALL  SINGl  (NC,N2*KSYl,KTcl,t.Tb2,CHJkDC, 

1  SwttPi*b.»tcP2,'>wt.tP*L)lHtDi,ulHtD2#OIHcO» 

2  ZS,XLE,YLE»XC,XZ»XZZ,YC»YZ#Y2Z, 

3  Z*Cl,C<:»C3i£i»t2,!.3#L-'f,E.^>lN0) 

CALL  SUPF   (ND,Nt,NCf NX»nZ# ISYr,KSYM,<lcl»KT£2,SCALf 

1  YA^jACZjZi^ACffCjLOPT.lkAlL^XSjYSfNP, 

2  IT<^l,ITc2»IV,bc,Zl,XP,  YP,D1#02,C3,X,Y,  IND) 
.0)  GL  Tu  291 


{ IND.Ft, 
TO  151 


IF 

GO 

2!>1  Kl 

K2 

QG  252  M«l 
WklTi-  {^) 
LG  2t2  K-1 
bUFFEK  IN 
IF  (LNIKN 

262  WRITE  (^) 
REWIND  Nl 
WfilT^  (<») 
END  FILE  ^ 

252  CONTINUE 
REWIND  ^ 
CALL  SSk»TC 
IF  (IbTOP. 
JIT 

!»-  (MT.LT 
GO  TO  161 

2«1    REWIND    ^ 
GO    TO    151 

291  WRITE     (IwR 
wmTE     (IwP 

292  F0R»1AT(24h 
GO    TC    1 

301    IF    (KPLQT. 

STOP 
500  FORMAT(IX) 
5iO  FoPIATCbFl 
530  F0R1AT(20A 
600  FORiATdHl 
61D  F0RiAT{F12 
620    F0R^AT(8fcl 


»    KTEl       -1 

•    KTe2      +lTfc2(KTE2)       -NX/2 

#3 

\XfK'(,KZ,>i^,Ki,KZ,l  IT 

(M,l)     (G(*,i#a)*G(hX#MY,  1)  ) 
1 ) .GT.C. )    GO    TG    2ol 

< (G(I»J»1)* I-l/MX), J»1,^Y) 

(F0(K),K-M,K2) 


H(l,iSTQP) 

EC!)    GO    TC    161 

»    0 

.KIT.ANO.AfcS(DG)  .GT.CUV.ANO.ASSiOG)  .LT.IO.)    GO    TO    I'll 


1T#60G) 
IT»292) 
Ofc<AO    LaTA»SPLInE    FAILURE) 

GT.O)    CALL    PLOT{0.»C.,99s») 


0.6) 

^) 

) 

.^#7F15.^) 

5,5) 
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630    FGRMAT(1HC»20A^) 
6^0    FCRiiATdb,  7Iii») 
6^C    FGRMATdX,  3^I<») 

660    F0RMAT(I10,E15.tJ»3I^,ti5.5,-3lAOF10.5,liJ) 
670    F0«i1AT(atl5.^>2Fl^.^) 

700  FOR^ATliSHOCGf-PLTING  TiMt  »  FIG  .i#  lOH    ScCCNDS) 
END 


StBROUTlNfc   GEOh* 


11 


1 

2 
GEQMET 
DIMtNS 

1 
IkEAD 
IwRIT 
KAD 
READ 
fcEAO 
IF  (FN 
KSYM 
NC 

WRITE 
2  FORMAT 

1 
WKITE 
WRITt 
SwEEPl 
S  w  F  t  P  2 
SWEEP 
DIHEDl 
DlHcC2 
DIHED 
ISYMC 
XlEO 
CHQRDC 
K 

READ 
READ 
ALPHA 
IF  (K, 
READ 
READ 
NU 
NL 
N 

kEAD 
kEAO 
READ 
UU  12 
PtAO 


PIC 
ION 
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DEF 
X 
S 

«  5 

■  6 

■  b 
(  IRtAD, 
(IREAD, 
C.L1.3. 

«  Z 
«  F 
(IWFlTf 
{15H0 

iSH 
(I.-RIT* 
(  IwRIT* 
S 

s 
s 

0 
D 
D 
1 
0 
0 

=  I 

(IREAD, 
( IREADf 

*  A 
GT.1.4N 
(  IREAD, 
(IRfcAC, 

-  F 

-  F 

•  N 
(iPEAD* 
( IRt AD, 
( IfiEAO> 
I»NL,N 
(IREAD, 


(Nb,.>(C,NP,ZS,XS,YS,XLt,YLE,SLOPT,TRAIL,AP,YF, 
SWtfcPl,SwttP2,SwEtP,0lHt01,0lHtD2,DlMtJ, 
XTtO,CHOKLO,ZTIP,lSYiO,KSYM) 
INITICN  OF  wING 

S(ND,i),Y3(N0,l),Zj(l)»XLE{i),YLE(l), 
LDPT(i),rKAIL(i),XP(l),YP(i),s?(l) 


7.2 

5  CO 
510 
)  k 
SY^1 

NC 

2) 

$ 

I 
61C 
61C 
WEE 
wEE 
WEE 
IHE 
iHt 
IHE 


9577v5130b23 
) 

)     ZSYM,FNC,SwLLPl,SwtEP2,i).>EEP,DHEri,DlHi:J2»tlHE0 
cTuPN 


NEEP(i)       ,IW  SWEEP(2)       ,15H 

iHtD(l)        ,lt?h  OiHzu(Z)        ,15H 

)    XL,  fL»CHUR[ , Thick,  *.L 

)     bw£cPl,SWtlP^,i«fcEP,CHE01,DlHcU2,LlHED 

Pl/KAD 

F2/RAJ 
P/RAD 

Dl/R/0 
U2/KA0 
P/RAD 


FINAL  i)»»EEP  , 
FINAL  UIHbD  ) 


5  00) 

bit)  ZS(K),XL»YL,CHQKu,THlCK, AL,FSEC 

L/RAD 

O.FitC.fci-,0,)  oL  TO  31 

5Ct  ) 

51C)  YSYh,FNU,FNL 

NU 

NL 

U   tNL   -I 

5  CO 

:?1C)  TRL,iLl,XSlN&,YSlNG 

500) 

51C)  XP{ 1) , YP(  I) 


58 


L  •    Hi       *1 

IF     (YSYM.GT.O.)     GZ     TG    ID 

Rt-AD       (IRtAD»t)CC) 

CO    l<i    I-1*NL 

PEAD       (IPfAD»i;IC)     \^t■l,OUl^^ 

J  •     L        -I 

XP(J)  "     VAL 

i^    YP( J)  -    DUK 

GO    rC    21 
li     J  =    L 

DL    11.    I-NL#N 

J  »    J       -1 

XP  (  J  )  ■     XP(i ) 

16    YP( J )  =    -YP( I) 

21    wkITt     (  Uf<iT,fcOG> 

WkITt     (UP1T,22)    ZS(r<) 
i?    FCRIATClf  nOPKCiFILt     Al     Z    »    *FiC.b/ 

1  itHO  IE    ASolE       tl3h  TE    SLOPE       »15H  X    ilNG 

2  ISh  Y    S]^G  ) 

wklTL     (UklT,fclC)     TFl»SLI»X:>1NG»  TMNG 

wRiri     {IwRIT»2<t) 
24    FGR^Al  (15H0  X  »  Ipr-  Y  ) 

UD    26     I-1,N 
26    Wkirf     (lWRIT»eU)     XP(1),Y0(I) 
3i     SCALF  =    ChLF  (?/(  XF  (  1)       -XP{NL)) 

XLE(K)  -    XL       ♦(XSIinG       -XP(NL  )  )''rHICK*sCALt 

YLF(^)      =  YL   -KYSING   -Y  P  (  NL  )  )  ♦  TH  ICK  ♦:)C  Al  E 

XX  "  XF{^L)   ••(XSi.MC   -XP(NL)  )*TH1CK 

YY  =  YP(M)   ♦(fbl.MG   -YP{M))*mCl< 

CA  ■  CaS(ALPHA) 

SA  «  SIN(ALPHA) 

LL  32  I"1*N 

XS(I»K)    =  SCALE*( (XP( I )   -XX)4CA   ♦ TH i CK ♦ ( YP ( I )   -YY)*6A) 
32  YS(I*K)     -  SCALt:*(ThiCK*(  YP(  I)   -YY)»CA   -(XPCl)   -XX)*SA) 

SLOPT(K)   '     IHICK»5LT   -T AN ( A LPHA ) 

TRAIL(K)   »  TmJCK*TkL/KAL 

NP(K)       =  N 

XTEO       -  AMAXKXTtOfXiC  i,K  )  ) 

CHQRDO     «  ArAXKCHLRUU^CHCiPD) 

IF  {YSYM,LE.O..Lk.ALF-iA..N  =  ,C.  )  iSYfO  «  0 

WRITE     (I».PIT,52)     7S{K) 
52    FCR«AT(27HOS-CTiON    utFi^ITIU^    AT    Z    =    /FIG. 5/ 

1  It^HO  XLE  ,ie>h  YLt  tl'jH  CHQKU 

2  13HTriICK^ci5    FATIC/UF-  ALPHA  ) 
UKlTf     (iWRlT,61G)     XL*YL*CHCRD*TH1CK,AL 

K  =    K       +1 

IF    (K.LE.NC)     GC    TO    11 
ZO  •    .5*(Zb(l)       ♦ZSINO) 

IF     (KSYM,Ne.3)     ZO    «    Zi)(i) 
DQ    62    K=1,NC 
62    ZS(K)  -    Zi(K)       -ZO 

ZTIP  «    ZS(NC) 

RETURN 
f>0C    FuR.-iATdX) 
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510    FC»tM*T(8F10.6) 

600    FURMATdHl) 

610    F0RMAT(Fli.A,7Fli).^) 


S03RCUTlNt       COCFD 


1 
2 


(KS 


SfcTS    U 
DirlbNS 

DX 
DY 
KY 
DZ 
ZC 
Kl 
K2 
IF 
DZ 
Zt 
Kl 
K2 
AX 
AY 
AZ 
BX 
bl 
XMAX 

zriAx 

SY 

SCAL 

i>CALZ 

V2 

Wl 

W2 

S73 

BBX 

ABX 

CBX 

A6BX 


DC 

DD 

B 

IF 

A 

AS 

C 

DC 

Dl 


12 


(AE 


P  STRFTCHt 
ION  ACCl 
Z(l) 
«  2./N 
«  i./t 
«  NY 
»     2./N 

■  1. 

m      1 

»    NZ 
YM.tQ.O)    G 

■  l./N 

•  '). 
=    3 

m     N|Z 

-  .;> 

«   *i 

•  0. 
=    0. 

■  .62t 
=    .625 

■  .  5 

«    XT£C 

-  ZTIF 
=  (DX/ 
»    SCAL 

«  (  Wl* 
«  SORT 
=    -tJX* 

-  1.    - 

«=     (  l**. 

•  4BX 
SORT 

I=2»NX 

»     (  I 
=    1. 

S(OD).GT.X 
"    C6X 
=    SORT 

■  ABXf 
=  ABX* 
»     4S/C 


(NX#NY#NZ»KSYf*XT-0,ZTIP*XnAX»ZMAX# 

SY»bCAL»SCALZ»AX>AY»AZ, 

AO»Al,A2,A3,bO»Bi,b<:,33fZ,Cl*C^,C3) 
PAKAtjOLK    AND     SPASxISc    CJOKClNATtS 
Al(l),A2(l)#A3(l)*dO(l)*bl(l),6?(l),&3(i), 
l(l)>C2(i),C3 (1) 


0 

)» 

fZ 

X 

Y 

♦  1 

Z 

-D 


0    TO    1 
Z 


♦  2 


/(  .eCCjl+KMAX  +  XMAX) 

/(i.f  O00ol*Z»'/iX  ) 

DY)**2 

/SCALZ 

QX/CZ)**^ 

(73.) 

S0«'T(3.*(7.    +    S73))/((l.    +    S73)  ♦Xf1AX**d) 

h<PX»SQRT(  (7.     +     n73)  /  12.  )*X1AX*»3 

+     S73)*XMAX*XMAX/12. 
♦     c*fiX»(3.*C«X    -     '(.♦XHAX^Xf  AX)*X1AX*XMAX/ 
(CBX    -    XriAX*XhAX) 

-1)*['X       -1. 

MAX)     GO    TQ    ]3 

-    OC*DO 

(A) 

AS    +    d9X*(3.*CBX    -    <t.*OD*CO)*DO*00 

DD    +    BBX»AS*DD*»3 
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13 


i<. 


12 


2Z 


02 

GO 

IF 

A 

C 

D 

CO 

Dl 

02 


AO(I  ) 

Aid) 

A2(I) 

A3(I) 

DC  (Ll 

DO 

A 

C 

D 

Di 

t)0(  J) 

bl(  J) 

bd(  J) 

e3(  J ) 

eez 

A6Z 

cez 
Aedz 


TO  1^ 
(DD.LT.O. 


J»2# 


>,2    K-2*K? 


(ABS(DD) 


33 


TO  34 
(DD.LT.O. 


34 


32 


Uu 

to 

8 

li- 

A 

AS 

C 

00 

Ll 

02 

GG 

IF 

A 

C 

u 
oc 

01 
02 


Z(K) 
Ci(K) 
C2(K) 
C3{K) 

RtTURN 


Q^X^CC-jH^C-b.+C  jX  +  Is'.+LO^yO)  -  12  .  ♦aC»*'» )  ♦DO/ (  A*C  ) 

)  8  -  -1. 

1.   -((OC   -B*XrAX)/(l.   -Xf^AX))»*2 

A*»AX 

(AX    +AX    -1. ) ♦(!,    -A) 

'3*X^AX  ♦  AfHX«(Or  -  P*X1AX)/C 

A*C/(  (  i.  ♦  0)*AfibX) 

-(AX   ♦AX)«(DL)   -R»X1AX) 

»(J.   ♦D)/((l.   ♦D)*A*(1,   -X^*AX)*♦2) 
DO 

.;)*(i/ox 

Dl  +  Dl  ■< 

. t*rx*D> 
(Kr  -j)»or 

i.   -00»lO 

4»*A  Y 

(AY   ♦AY   -l.)*(l.   -A) 

A*C/( (i.   +J)*SY) 

iY*DD/C 

.t;*Ll/DY 

1 1  ♦  r  I  ♦  V  2 

-AYvCn»L)T«(  3.   ♦C)/((l.   ♦DJ^A) 

-e2*SUKT(3.*(7.  ♦  S73))/((l.  +  S73)*Z1AX*»3) 

1.  -  bijZ*SwRT(  (  7.  ♦  S7j)/lt.  )*ZMAX**3 

(!«;.  ♦  S7i)*Z1AX»ZVAX/12. 

AbZ  ♦  bdZ*(3.»CBZ  -  4.*ZrAX*Z«AX)*ZnAX*ZMAX/ 

SwRKCbZ  -  ZnAX»ZMAX) 

(K       -KlJ*uZ       -ZO 
1. 

■  GT.Zr.AX)    GL    TQ    i3 
CtZ    -    UD*0O 
SC<iT(A) 

ABZ*AS    +    BaZ*(3.*CBZ    -    4. ♦DD+DD ) ♦DD*OD 
AdZ*OC    ♦    B8Z*AS*rL**3 
AS/C 
BtfZ*(C«Z*(-6.*CBZ    +    19.»DD*0D)    -    12. ♦0D**4  )  *Du/ (A^C ) 

)     «    •    -i. 

i.      -((OD      -t3*ZM*X) /(I.       -Z«AX))**2 

A  +  +  AZ 

(AZ       ♦AZ       -l.)*(l.       -A) 

8*ZfAX    ♦    A!JttZ»(OL     -    H*ZrtAX)/C 

A*C/( (1.    +    D)*ABeZ) 

-(AZ       ♦AZ)*(uO      -b*Z.1AX) 

♦  (j.       ♦D)/((i.       ♦OJ^A^d.       -ZHAX)*»i:) 
SCALZ*00 
.5*L1*W1/0Z 
Dl+Dl*w£ 
.;?*DZ*02 


61 


SL'BRCUTiNt 


1 
2 

3 


GtNfcRATES 
CIMESSIQN 


1 
? 

bii    2  K« 
t^(K) 
2  65(K) 
Kl 

IF  (KSYM.E& 

M 

K2 

KTtl 

11  DO  12  K«K1, 
IF  {Z(K).LT 
IF  (Z(K)  .L6 

12  CGNTINUE 
B 

Si 

S2 

Tl 

T2 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

CALL 

S 

Si 

T 

Tl 

T2 

XC(2) 

YC(2) 

IF  (KSYM.N 

N 

00  22  K=K1 

11 

A 

XC(K) 

YC{K) 

XZ(K) 

YZ(K) 

XZZIK) 
22  YZZ(K) 
31  N 


0. 
2 

NZ 
.0) 

3 

NZ 

3 
K2 
.ZS 
.ZS 


GO  TG  11 
+  2 


(1  )  )  Klci  ' 
(NO)  kT£2 


♦  1 


SPLIF 
INTPL 
INTPL 
INTPL 
SPLIF 
I^TPL 
INTPL 
INTPL 


CHOPOO 

TAN( SWct Pi) 

TANiS-^Ef  P2) 

TAN(DIHtCi) 

TAN(DlHtD2) 
(i*NC>ZS*XLif.:i»t2,F3»i»Si^l*S2»0»0.»iND) 
(KTEl»KTEi*Z*XC»l*NC,Zb»XLE»cl#b2»E3»0) 
(KTEl*KTc2*Z»XZ*l>NC»ZS*cl#£2»t3*t4»0) 
{!<>TLl»KT£2»Z*XZZ»:*NC/ZS»t2»i3»LA»t5»0) 
(l»NC»ZS»YL£»fcl>E?,E3»i*Tl»l»T<:fU#0.*lN0) 
(KTei*KTE^,2*YC»l*NC,ZS*YL£»El,h2#t3,0) 
(KTEl>KTE2fZ»YZ»l#NC*Zi>*tl/E2»E3»cA,0) 
(KTEl*KTc2»Z»YZZ»l*NC»ZS»f2*t3»b<t»t5»C) 

B*TAN(S».EtP) 

&*S1 

b*T*^N(OIhED) 

6*T1 

B  +  T2 

3.*(XC(3)   -XC(H))   +XC(t') 

3.«(YC(3)  'fZi'*))       +YC(t/) 
,0)  to  TG  31 

KTFl   -1 
N 

{Z(K)   -2(KTtl))/e 

fcXP(ZZ) 


XC(KTtl) 

YC(KTEl) 

(S   +(bl 

(T   +(T1 

(SI   -S)*A/(e*6) 

(71   -T)*A/(b»e) 

KTE2   +1 


♦S*ZZ   -(SI 
+T»ZZ   -(Tl 
-S)*A)/h 
-T)*A) /B 


-S)Mi.   -A) 
-T)*(l.   -A) 
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32 


Ll    32    K.N,K£ 
ZZ 

A 

XC(K) 

YC(K) 

XHK) 

Y2(K) 

XZZ(K) 

YZZ(K  ) 

PtTUf-M 


( Z(K )       -2(K  rt2) )/B 
EXP(-ZZ) 


XC(KTEi) 
YC(KTL2) 
( b  +(S2 
{ T       +(T2 


♦S*ZZ       +(S2 
♦  MZZ       4(12 
-.S)*A)/E 
-T)*A)/E« 


-(SZ       -S)»A/(B*B) 


-T)*(l. 


-A) 
-A) 


StBRPUTlNt   SbPF 


1 
2 


1ST 
INT 


ckPQL 
dNSID 


1 
2 
3 


21 


23 
2^ 


32 


PI 
TYA 

Si 

DX 

LX 

MX 

MZ 

IVO 

IVl 

Du 

HE 

ITE 

DG 

IV( 

SC( 

K 

K2 

K2 

M 

6  2 

IF 

fc2 

Rl 

C 

cc 

OG 

IF 

IF 

CON 

lit 


ATCS 
ATI'oN 

N  li 
X 
X 

I 
3 

T 

• 
2 
N 
N 
N 
1 


2  (<»1, 

i(K) 

2(K) 

2  I"l/ 

I*K) 

I»K) 


(ZS(K2) 


.■"-APP 

1^ 
0(NE 
C(l) 
P(l  ) 
V(Nt 
.Mi 
An(  Y 

.  /KA 
X/2 
X   + 

Z   ♦ 
-I 

1   - 

X 
X 

2 

Thi 


(  N 
Y 
I 
Lb 
LlN 
»1) 
»YC 
f  YP 
>1) 
5'?2 
A>.) 
Al 


C#NE*NC*KX*NZ»ISY,^,KSYH,M£l,KTe2*SCAL» 

A*»AO,Z»ZS*XC,YC»SLbPT,TKAIL»XS#YSfNP* 

lfci.fiTt2»lV,SC»Z0,XP»YK,Ul,D2,  0  3*X,Y,ihU) 

WlN(j  SORFACt  AT  MEiH  PulNTS 

t A^  IN  PHYSICAL  PLANE 

*Ai(ND»i)*YS(NJ,i)»Zb(i)*SLJPr(l)*TRAiL(l)# 

(  1)»A0(1),Z{1)»ZC  (1  ).X(1),Y(  ]), 

(  1)  »DI(1)  M  2(1  )  ,LH  I)  t 

»NP(i)»lTtl{i)/ITE2(i) 

6535t97y 


♦  1 
1 
3 

SY^ 
ISYh 


-ISYfl   -ISYr 


32  I- 
(  (A0( 
(  (A0( 

TINUE 
i(K) 


2» 
I) 
I) 


NX 


2   ♦! 

2   -1 

• 

-Z(K))  21#2'j*23 

{Z(K)   -ZS(Ki) )/(ZS(K2)   -ZS(Kl)) 
i.   -R? 

R1*XS(1»M)   ♦R2»XS(1»K2) 
SQPT{(C   +C)/SCAL) 

.r.*DX).LT.-CC)  11  »  I  *l 
•.5*0X).LT.CC)  12  «  I 


II 
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<.l 


*E 


nt2(K 
cc 

ZO(K) 

KK 

P 

N 

U 

00  (^2 

X(I) 

ANGL 

U 

V 

DO  4^ 

IF  (i«, 

ANGL 


I«2*NX 


EC.C 


=  S 
.  ) 


j(U)/CC 

(K)   -lYAw*(XC(K)   +Si«AO( IZ)*A0(I2)) 

1 

1 

P(KK  ) 

U«T( X5( 1>KK)/C)/CC 

♦  A0(  I) 
1   +P1 


43 


44 


U  •  X 

V  =  Y 
R  -  5 
XP( I)  »  R 
YP(I)  »  ft 
GC  TO  44 

ANGL  =  PI 
L  -  -1. 

V  «  0. 
XPd)  -  C. 
YP(I)  =»  C. 
CONTINUE 

ANGL 

ANGLl      « 

ANGLE 

ANGLl 

ANGL2 

Ti 

Ti: 

CALL  SPLIF  ( 

CALL  INTPL  ( 

XI 

A 

B 

ANGL 

U 

V 

M 

UO  52  I  =  2*M 

XX 

D 

YY 

R 

ANGL 

b 

V 

R 


Of'K  XS(  I*KK)**2   ♦YS(i»^K)*»2) 
GO  TO  4: 

NGL   ♦AI AN2( (u*Y:(  I»KK) 
(U»Xb(  I*KK) 
S(1,KK) 
b(l»<K  ) 

QRT((t<   +R)/jCaL) 
♦COS{ .S^ANGL) 
*S1N(  ,5*ANGl) 


•»V*YS(  i,KK))  ) 


4T<sh  (  SL(JPT(KI<  )  ) 

aTan( YS(l»^^)/XS(x»KK) ) 

ATAN( Yb(^#^K) /Xi(N*KK) ) 

ANGL   -.5'*(AN&Li   -TPAIL(KK)) 

ANGL   -.i*(ANGL2   +rf<AlL(KK)) 

TAN(ANGlI) 

TAN( ANtL2) 

l,N,XP,YP,L)i*O2»L3>i»Tl»i,T2,C*C.>l.N0) 

Il»I2*X#Y»i»rfXP,YP,[;l,J^»C3»C) 

,25*X5(  1#K(<) 

SLCPT(KK)*(XS(l»t<K)   -XI) 

l./(XS(  1»KK.)    -XI) 

PI    +P1 

1. 

0. 

II   -1 

.5*SCAL*X( I)**2 

•J*(XX   -XI) 

Y5-(1»KK)   ♦A*ALU{  (L  )/0 

SQR1(XX**2   +YY*»2) 

ANGL   +aTAN2((U*YY   -V*XX)f{U*XX   +V*YY)) 

X  X 

YY 

SGR1((K   +k)/$CAL) 
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SLLPrCht^  )  ♦(XS(N»f  K)   -XI) 

i./(XS(N,KK)   -XI) 

J. 

1. 

0. 


.'♦SCAL*X( I)»*^ 

ft* (XX  -XI) 

YS(N,KK)   ♦A*ALaC-(P)/D 
SQPT(XX**t   +YY**2) 

ANC-L   ♦AIanZC  (0*YY   -V+XXJfCJ+XX   ♦V^YY)) 
XX 
>  Y 

SOPT((P   ♦RJ/SCoL) 
K»SIN(  .'^♦ANGL  ) 
P*C»CC*CC 
X 
SO(I»K)   ♦Q*Y(I) 

2)  Gc  TO  n 

K2 


5?  Yd) 
A 

d 

ANGL 

U 

V 

r 

XX 

b 

YY 

R 

ANGL 

U 

V 

P 
5A  Yd) 

C 

00  a?    I»2» 
62  SC(1»K) 

IF  (KK.EO.K 

KK 

P 

Gu  TC  -tl 

71  LO  72  I'll, 12 

72  IVd»K)    »  2 

r  -  il   -I 

D(J  7<.  I«2»l^ 

Z2  •  Z(K)   -TYA**(XC(K)   *il«AOd)*AO(I)  ) 

If  (ZZ.Gt .ZO(KTPl) )  iv(I>^)  «  IVJ 
7't  CCNTINOF 

PI  «  12   ♦! 

DQ  7t  I«>"#NX 

Z2  ■  ZC^)   -TYA>,«(  XC(K)   ♦SI*AOd)*AO(l)  ) 

IF  (ZZ. Gh. ZO(KTtl)  )  IVd»K)  «  IVw 
76  CCNTlNUf: 

K2  ■  K2   -1 

K  ■  K   •»! 

IF  (K.Lfc.KT62)  GO  TO  21 

Kl         -  2 

K2         «  NZ 

IF  (KSYM.iQ.O)  GO  TC  61 

Kl  -  ^ 

K2  »  NZ   <-2 

81  00  32  I-2»NX 

ZZ         «  Z(K)   -TYA»<*{XC(K)   ♦Sl»A0d)*A0d)) 

IF  (Z2,LL.ZS(NC).ANL;.ZZ.Gi.ZC(KTFi)  )  l\rf(I,K)  -  IVO 
62  CONTINUE 

K  =  K   +1 

IF  (K.LL.K2)  GC  TO  fl 

N  •  KTb2 

IF  (YAW.LE.O. )  GO  TO  93 

10  ■  ITcl(KTF2)   *i 

DU  92  I«IC*LX 
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N 

ZO(N) 

I 

ZC(KTfcl-l) 

Z0(N-»1) 
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(  .5*(X(1*1) 
.&*( Y( I+i) 
.  !J*(CP(  I  +  l) 
-CPA+DX 
CPA*uY 
CL   +DCL 
CO   ♦DCO 
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=  CL^CUSC  ALt^hA)        -CD*  SlN(  ALPHA) 

■  CL^SINiALf'HA)        4CL*C0S(ALPmA) 
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CALCLLA 
OntNSI 

CL 


TPS 

ON 


TO 
r 

Z 

C 
G 
0 

u 

G 
0 

K 


CTFU^.{^Ui»KT•:^#CHLl^CfJULoCD#Sv:M*Z*XC» 
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Z(35)»C1(33),C2(35),C3(35)» 
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MY         '    NY   +2 
MZ         »  NZ   +3 
MXO         =  NX/2  *l 
MYO         •  NY/2   +2 
KZO         »  NZ/2   +1 
K  ■  i. 
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»  K 

112  J"i>MY 

112  1=1, MX 
G(I» J,2)   ■  .5*(G( I, J*l) 
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»K  +  1)       -SOd>K-l) 

)*Cil 

)*DSK 

C(i»I\/(I,K)) 

-K        ♦AJ(l)*AGd) 


+  SCd,K)*SO(I#K) 


D^XZCK)       -i(  d,K)*YZ(K) 

I)*Y1(K)        +SUd*K)»XZ(K) 

X       -eZ       ♦Fh+SZ 

♦  SX»i)X       •«-H*HZ*HZ 

+H^*Z*HZ 

1#KY,2)       -G(I-1#KY,2) 

KY»3)       -Gd.KY^l) 
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This  report  was  prepared  as  an  account  of 
Government  sponsored  work.   Neither  the 
United  States,  nor  the  Administration, 
nor  any  person  acting  on  behalf  of  the 
Administration: 

A.  Makes  any  warranty  or  representation, 
express  or  Implied,  with  respect  to  the 
accuracy,  completeness,  or  usefulness  of 
the  information  contained  in  this  report, 
or  that  the  use  of  any  Information, 
apparatus,  method,  or  process  disclosed 
in  this  report  may  not  Infringe  privately 
owned  rights;  or 

B.  Assumes  any  liabilities  with  respect  to 
the  use  of,  or  for  damages  resulting  from 
the  use  of  any  information,  apparatus, 
method,  or  process  disclosed  in  this 
report . 

As  used  in  the  above,  "person  acting  on  behalf 
of  the  Administration"  Includes  any  employee 
or  contractor  of  the  Administration,  or 
employee  of  such  contractor,  to  the  extent 
that  such  employee  or  contractor  of  the 
Administration,  or  employee  of  such  contractor 
prepares,  disseminates,  or  provides  access  to, 
any  information  pursuant  to  his  emplojnnent  or 
contract  with  the  Administration,  or  his 
employment  with  such  contractor. 
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